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TECHNICAL PAPER 


A NUMERICAL SOLUTION OF THE NAV I ER -STOKES EQUATIONS 
FOR CHEMICALLY NONEQUILIBRIUM, MERGED STAGNATION 
SHOCK LAYERS ON SPHERES AND TWO-DIMENSIONAL 

CYLINDERS IN AIR 

SUMMARY 


The complete Navier -Stokes equations are solved along the stagnation 
streamline in merged stagnation shock layers on spheres and two-dimensional 
cylinders using an iterative finite-difference numerical procedure known as the 
accelerated successive replacement method. The fluid medium is chemically 
reacting air consisting of seven species. Velocity components, thermodynamic 
properties, species mass fractions, and wall heat transfer rates are computed. 
This report is intended as an explanation of the method and as a user’s manual 
for the computer program. 


I. INTRODUCTION 


An aerospace vehicle ascending or descending through the Earth’s 
atmosphere traverses several flow regimes from the continuum boundary layer 
regime at low altitudes, through the transitional regime at intermediate alti- 
tudes, to the free molecular regime at very high altitudes. The character of 
the flow field changes drastically from the boundary layer regime to the free 
molecular regime, and no single computational approach is valid throughout 
this range. The broad transitional regime maybe divided into several sub- 
regimes as suggested by Hayes and Probstein [1]. Consider the typical 
spherical nose of an aerospace vehicle. In the boundary layer regime viscous 
effects are primarily confined to a thin boundary layer, the bow shock can be 
treated as a discontinuity, and a region of inviscid flow exists between the shock 
and the boundary layer. However with increasing altitude, the shock wave and 
boundary layer thicken and eventually merge into a single viscous layer called 
the shock layer. The flow regime in which this occurs is called the fully 
merged shock layer regime which is the condition treated in this report. 



At the great speed that a vehicle reenters the atmosphere, the tempera- 
ture near the body becomes extremely high, especially in the stagnation region. 
Therefore, the air in the shock layer dissociates and ionizes. For an accurate 
description of the flow field, one must account for these real gas effects. Also 
an accurate estimate of the ionization level is needed for radio communication 
purposes. Therefore, a flow field model including finite rate chemistry is 
required. 

This report describes a method for computing flow properties along the 
stagnation streamlines of a sphere and a circular cylinder transverse to the 
flow. Heat transfer rates are computed at the body surface. Although this 
method is limited to the stagnation region, it still provides valuable design 
information because maximum heating rates usually occur at the nose. The 
computational method was developed by Jain and Adimurthy [2] for an ideal gas. 
The method uses the full Navier-Stokes equations to describe the flow in the 
entire shock layer from the surface to the freestream. The boundary conditions 
at the wall are provided by slip velocity and temperature jump equations. Using 
the concept of local similarity, the governing equations are reduced to a system 
of nonlinear, coupled ordinary differential equations. Numerical solutions are 
obtained for points on the stagnation streamline using an iterative finite - 
difference procedure known as the accelerated successive replacement method. 
The applicability of this approach and the failure of thin-layer theories for the 
merged shock layer regime is discussed in Reference 2. Nonequilibrium 
chemical reactions were included in this method by Kumar and Jain [3] using an 
air model with seven species and six reactions. Hendricks [4] developed 
surface slip velocity and temperature jump equations for a multi -component 
gas, including the effects of wall catalysis, to use with this model. Additional 
modifications have been made in this report, principally by including the two- 
dimensional cylindrical geometry and using a multi-component gas model to 
compute viscosity. 


II. ANALYSIS 

A. Formulation of the Problem 

1. Approach . In the present analysis the full Navier-Stokes equations, 
with nonequilibrium chemistry, are solved through the merged stagnation shock 
layer from the freestream to the body. The slip conditions at the gas-wall 
interface include the effect of wall catalysis and a multicomponent, nonequilibrium 
gas flow. 
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The thin shock layer assumption is not made in the present analysis. The 
full Navier-Stokes equations with chemically-reacting, nonequilibrium air are 
solved through the merged shock and boundary layer. This allows the shock 
wave to develop within the computational domain. A seven species air model 
is used. The species considered are N 2 , 0 2 , NO, N, O, NO + , and e - . For air 
dissociation and ionization, the rate expressions recommended by Wray [5] are 
adopted. Prandtl number, Pr, and Lewis number, Le.., are taken to be 0.75 

and 1.4, respectively, for the cases computed in this report. The viscosity of 
dissociated and ionized air is obtained from a simple summation formula for a 
mixture of hard spherical molecules; using Hansen’s collision cross sections [6], 

Solutions are obtained by using the local similarity concept to reduce the 
governing equations to a set of nonlinear, coupled, ordinary differential equa- 
tions. This set of equations is integrated using a finite difference method known 
as the accelerated successive replacement method. It is important to note that 
the first order local similarity assumption is a good approximation near the 
stagnation streamline, at least for Re^ 10 [7], Reference 7 gives a 

thorough discussion of local similarity. 

2. Governing Equations . The nonlinear, coupled ordinary differential 
equations governing the flow of a multicomponent gas near the stagnation stream- 
lines of spheres and two-dimensional cylinders are presented. The coordinate 
system employed is shown in Figure la. 

a. Basic Equations. The basic conservation equations and the ideal 
gas equation of state for a multicomponent, reacting gas mixture are as 
follows [8]: 

Global Continuity (of all species): 

— ^ + V • (p v) = 0 , (1) 

at 


Species Continuity: 

D Y i - 

P— =- = w - V • (p Y V ) i = 1, ..., NS , (2) 

Dt 


where NS = number of species in mixture (no summation on repeated indices). 
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Momentum: 


- D v 

P 

Dt 


NS _ 

- V p V . T + p Yj Y i\ 

i=l 


( 3 ) 


The second order tensor, f , is the viscous stress tensor and f is the body 
force per imit mass of species i . 

Enthalpy: 


-Dh 

P — 

Dt 


+ $ - V • q 

Dt 


( 4 ) 


The quantity $ is the viscous dissipation function and q is the heat flux vector. 
State: 


NS Y. 

P - pfi T ^ 


. . W. 

1=1 l 


( 5 ) 


b. Nondime nsional Equations. The basic equations are put in non- 
dimensional form by introducing dimensionless variables as follows: 


u = u/V^ 
v = v/V ro , 

P = P/P ro 
T = T/T 0ro , 


5 



P = p/p m v 2 

r C0 00 

h = h/V 2 

00 


M = M/M, 0ro * 

r = r/f b , 

*i - V (5 » 


where is the radius of the body, T 0co is the freestream stagnation tempera- 
ture, and jil^ is the coefficient of absolute viscosity evaluated at T . The 
nondimens ional similarity parameters 

Re„ = p V r. /jl n , 

000 K C0 00 b ^Oco 

Pr = C p/k , 

P 

Sc.. = p/p D.. , 

i] i] 

and 

Le.. = k/p C D.. 
i] P ij 


are also introduced. The basic equations are then simplified by assuming steady 
flow and Newtonian fluid, neglecting body forces, viscous diffusion stresses, 
thermal radiation, and thermal diffusion, and using Fick's Law of Diffusion. 

The equations are presented in cylindrical and spherical coordinates below. 
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Cylindrical Coordinates: 
Global Continuity: 


A <P rv > + ^ (p») = 0 


(«) 


Species Continuity: 


1 £ ♦ H £) - *. + -i- i U(^ ♦ A (JL. S S) 

H \ dr r be) i Re^ r |_9r\Sc dr) be \Sc r be) 


i = 1, . . . , NS 


( 7 ) 


Transverse Momentum (0 Direction): 


/ 9_u u 9u vu \ _ _ 1_ 9 P + 1 

9r r b9 r / r 90 Re 


Oro 


1 _JL 

r 5 " 0r 


„2 

r^ 


9 / u \ + 1 9v 
9r V r / r 90 


19/ L/l 9u v\ 2 

+ r 90 I \r 90 r/ " 3 

Radial Momentum (r Direction): 


19 , , 1 9u 

— — ( rv) + — — — 
r 9r ’ r be 


( 8 ) 


p f v 9v + u 9v _u^\ = _ 8£ + _JL_ |"l ( L 8v _ £ [l 8(rv) 1 

\ 9r r 90 r / 9r Re 0co r 9r \ I 8r 3 r 9r r 


9u 

90 


r 90 



p -] 


/ 1 


8_/u\ 1 9v 


1 / 1 

4 

r 9r \ r/ r be 


“ 7 V 


“ 


\ l 


2 

3 


19/ s 1 9u 

— — — ( rv) + 

r 9r r 90 


O) 
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Enthalpy: 


/ V ah + U8h\ , V 8P + U 8P + ^r l/M* + /l 8u + vy) 
\ 0r r 80 / 8r r. 80 Re Qco [\8r/ \r 80 r ) j 


/l 8 v [ 8u uV 
80 8r r J 


2 _J£ 

3 Re 


0® 


8 v | 1 8u t v 
8r r 80 r 


r Pr Re 


i — \j.Lsk) + ^(iL$k) 

Re 0 oo| 8r \ 8r / 90 \r 80/ 



NS 8 Y. 


NS juh. 

8 Y." 


+ (Le - 1) — 

£ 

+ - h ig 

l 

00 


1=1 


1=1 

_ 

i 


(10) 


State: 


. « T o® N v s Y i 

P P T r, 2 Aw. 

V ^ 1=1 l 

00 


( 11 ) 


Spherical Coordinates: 

Global Continuity: 


? £ (f r2v ) + m (pu sln 9 > = 0 • 


( 12 ) 
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Species Continuity: 


P 


/ 9Y. 9Y. 

j 1 u 1 

\ 8r r 80 


) 


w. + 

l 



1 

r 


1 8 


r 8r 



sxn0 8 0 



i = 1, . . . , NS ( 13) 


Transverse Momentum (0 Direction): 



Eadial Momentum (r Direction): 
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Enthalpy: 


/ 9h_ + u 9h\ = 9P | u 8P t n J/8v\ 2 [ /l 9u + v V 

9r r 36 J V 9r r 90 Re j\ 9r/ \r 9 0 r/ 




cot 6 \ '' 


V 


[liz tI iM 

I r 90 9r \r/ 


(A 2 

.2 f» 

W 

3 Re_ 


o® 


9 v 1 9u 
9r r 90 


2v , u cot 0 

2 . i ] 

(_L JL , 

( r 2 

+ 1 8 / 

( \i sin0 9h\ 

r r 

Ee 0 co Pr 1 

[ r 2 9 r 


r sin0 8 0 \ 

V r df >) 


NS 


9 Y.\ 

+ (Le - x) s SF £ h i 77/ 


, „ NS /uh. sin 0 9Y 

MLe-1)^^ EH 


r sin0 9 0 H. \ r 
i=l 


B Ji) 

90 / 


(16) 


State: 


(RT NS Y. 

„ ^ 0® v 1 

P “ _ LJ w 

V 2 i=l W i 

00 


(17) 


In the previous equations, Le„ is assumed to be the same for all species 

pairs i,j and is, therefore, replaced with Le . Similarly, Sc., is replaced 
with Sc. X '' 


The local similarity approximation for the zone near the axis of symmetry 
is given by equations (18) through (24). The validity of this approximation was 
demonstrated by Kao [ 7] . 
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( 18 ) 


u(r, 0) = u^r) sin 0 
v(r,0) = v 1 (r) cos 0 (19) 

h(r,0) = hj(r) (20) 

P(r,0) = Pt(r) (21) 

P(r,0) = P t (r) + P 2 (r) sin 2 0 (22) 

P(r,0) = Pi(r) (23) 

Y i( r »0) = Y 1 .(r) . (24) 


Equations (18) through (24) are used to reduce the governing equations to a set 
of nonlinear, coupled, ordinary differential equations which can be solved 
rapidly. 

A transformation is now made in the normal coordinate by defining 


00 00 


where r^ is the nondimens ional distance from the origin to the freestream and 
n^ is the nondimensional distance from the body to the freestream (Fig. lb). 
The values of r^ and n m are unknown a priori. They are determined as part 
of the solution. This transformation, equation (25) , keeps the body at tj = 0 
and the freestream at n = 1 . 

By substituting equations (18) through (25) into the governing equations 
and equating the coefficients of like functions of 0 , one obtains the following 
system of ordinary differential equations: 
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111 mill 


II II i mi III 


■I Min II null 


Cylindrical Coordinates: 


p i v i 1 + ^ 


v. 

t 1 


R wi r i i u i + u i m , 2 p 2 Re o~ 


1 + T ? n co / P 1 ( 1 + r l n J 


U l + v i K 7 1 U l/ 2 , ^ 

+ (l + ^ n co) +3 ( 1 + nn a) ) "nll+r^ 


3n ( 1 + nn ) 

m ' 1 m ' 


3 E V , X 

(P T + n v vM I — + 

4 Vi K 1 1 1 1 ; n co Ui n co 1+ ^ n a 


VV.K 7 

1 + T/n ro 2 f i i n co + 4(l + t ? n oo ) 


P 2 P 1 P 1 U 1 

— = + --■■■■■ (u, + v ) 

n « n co 1 + ^eo 1 1 



p,vjl’ 

-A aa ( 1+ ’1”J 2 - ( 1 + »”J 2 


00 


V 1 P 1 ^ v 


t2 

1 1 


n 


CO 


Re„ n z 

Qco oo 


2m 1 

+ (u, + v )‘ 

Re Am' 1 1 

Qco 


2 


3 Re 


Qco 


( 1 + r?n ) 

v + ( u + v } 

n l K 1 V 


( 1 + ^ cq ) 

Re 0 co Prn co 


( 1 + vnj 

K h i + “ — iw f'i + 


+ ( Le - 1) 


+ 


4^(1 +inj \NS 


n 


\ns 

)& h ‘ Y “ 


n(i +v njm 

A (h!Y' + h.Y") 
' l li l lr 


co i=l 


(30) 


p v Y T 
P 1 1 li 

n 


= w. + 


i (1 + tjnJ Re 0o . Sc 




n 


n 2 

CO 




n 2 

CO 


i = 1, . . . , NS 


(31) 



(32) 


Here, a prime denotes differentiation with respect to 17 . 
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Spherical Coordinates: 


P_1 _ v _[ 2n co 
Pi " " V 1 " 1 + 



u" Re o /v u T u 2 + u v \ 

II = 0” 1 ( 11 + 1 11 

n2 u \ n 1 + Tin / 

CO ^1 \ CO ' CO / 


2P 2 R V 

Pi ( 1 + 7 7 n c0 ) 


U 1 + V 1 
+ ( 1 + *?0 


Pi 


1 00 


3(1 + rjn m ) 


S. / 2 

" n co i 1+r ? n co 


J^a 

Pi n , 


v i 

A ( 1 + *? n J 


v i’ 3 Re o« V 1/ P I 

rf ' 4 777 ( V p iVi> fe + T 

00 CO ^ 1 CO\^lCO 


+ nn 


2n 


(U 1 + V l> 

( 1 + i? n J 


P 


PA 


2 (1 + t? 11 ,*,) 


n 


p i p i u i 

— + ■■■ 1 (u, + v ) 

1 + r ? n m 1 1 


(33) 

) 

(34) 



( X + ? 7 n oo) 

(35) 

(36) 
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p 1 v 1 h 1 

-1^ (l + -inj 2 = (1 + nnJ 2 


t 2 


v i p i *1 v i 


oo 


Re oo> n <£ 


4 ix. 


Re 


Qco 


(« x + \) 2 


2 p l 


3 Re 


Qco 


d + T ? n J 


n 


co 


v' ± + 2 (u x + v i ) 


U +T J n ra ) 

Re 0» Pr n co 


(l + 7?n ® ) , 

*i h i + — K h i + ' , i h ? 


+ (Le - 1) 


' H + V*J\ *3 
2fx H — J X h.Y* 

\ 1 n / l li 

' co 1=1 


+ 


#*1 ( X + ’ I 11 ®) ^ 

— 7 X (hrr* +h.Y") 

n2 .H, i li i lr 


co 1=1 


(37) 


P l V l Y ii 


n 


= w + 


i (1 + rjnJ Re Qco Sc 


2^i Y * fx Y’ (1 +rjn ) 
r l li *1 li v ' 


n 


n2 

co 


V 1+ ’!»J Y i\ 


n2 

CO 


i= 1, ...» NS . (38) 


<RT. NS / Y ' 

*! - ~l(^ 

V 2 1=1 \ W i/ 


(39) 
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Equations (26) through (32) and (33) through (39) constitute two sets of non- 
linear ordinary differential equations for cylindrical and spherical stagnation 
regions, respectively. Each of these sets contain 6 + NS equations and 8 + 2NS 
unknowns (u , v , P , P , p , T , h , a, w., and Y ). The required addi- 

tional information and equations are given below for the mass production rates, 
w., enthalphy, h^, and viscosity, p^, of a dissociated and ionized air mixture. 

These equations, like the equation of state above, are independent of the 
coordinate system. 

c. Air Chemistry. The freestream air is assumed to consist of N 2 
and O z molecules only. The air model used for the shock layer consists of 
seven species and the seven chemical reactions as follows [5]: 


0 0 + M + 5. 1 eV p- O + O + M 
2 k 


fl 

S rl 


f2 

N +M+9. 8eV — N + N + M 
k r2 


f3 

NO + M + 6. 5eV - N + O + M 
k r3 


f4 

NO + O + 1. 4 e V ^ O q + N 

k . 2 

r4 


f5 

N. + O + 3. 3 eV — NO + N 

k , 
r5 


f6 

N + O + 1. 9 eV ^ NO + NO 
22 k 

r6 
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f 7 


+ 


N + O + 2. 8 eV — NCT + e 


r7 


These reactions are written so that the forward reactions are endothermic; the 
net amounts of energy required to produce the reactions are given on the left 
side of the equations. The first three reactions are the neutral-particle 
dissociation-recombination reactions in which the energy of dissociation in the 
forward reaction is taken primarily from kinetic energy by means of a collision 
with a "catalytic” molecule M ; the chemical energy released in the recombina- 
tion (reverse reaction) is converted primarily to kinetic energy in a three-body 
collision involving a catalytic molecule M . The M molecule can be any of the 
six molecular or atomic species present in the air mixture. The quantities k^. 


and k . are the temperature -dependent forward and reverse reaction rates, 

respectively, for the jth reaction. The experimentally determined reaction 
rates recommended by Wray [5] are used, together with the concentrations of 
all constituents, to obtain the mass production rates, w., of each constituent [9] . 
The reader is referred to Reference 10 for details. 


d. Enthalpy. The specific enthalpy of a mixture of gases is given in 
dimensional form as follows: 


Ns 

h = P v + £ Y.e. , (40) 

i=l 

where 


e. = e°. + e^ + +e (41) 

l l T. R. v. 

ill 


and 


e? = specific energy of formation of species i at the reference 
temperature (zero absolute) 

e^, = specific energy of random translation 
i 
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e = specific energy of rotation (for diatomic molecule) 

R. 

1 

e = specific energy of vibration (for diatomic molecules) . 
i 


We assume that the energy of electronic excited states is negligible, that 
the rotational State is fully excited, but that the vibrational state is partially 
excited. Values for these quantities are given by Vincenti and Kruger [9]: 


S T. = I * 
1 


(42) 


V = R i T 

x 


(43) 


e v " © /T 

1 Vi 1 
e - 1 


R. © 

1 v. 

l 


(44) 


where 


© = characteristic temperature for vibration of diatomic species i 

v. 
i 


© = 2270 K 


© = 3390 K 


© = 2740 K 

v 

NO 


© = 2740 K . 

V NO + 


18 



Equation (40) now becomes 


c NS Y. Y. , 

5 = o «T E w E £ T 


© 


i=l 


W. 

1 


. W. 

l l 




v. 

X 


NS 


(for 
diatomic 
molecules) 


© /T 

v. 

e 1 - 1. 


2 *! • 


i=l 


(45) 


Nondimensionalizing equation (45) , we obtain 


h = - 


_ RT T 
5 0“ 


2 - o 

00 


NS Y. 

E — 

4 1 w. 

1=1 l 


+ 


, NS 

E Y. e^ 

x"r 2 , . 1 1 


V' i=l 

CO 


(RT 


Qco 



(for 

diatomic 

molecules) 


(46) 


e. Viscosity. The Sutherland formula for the viscosity of air gives 
acceptable results at moderate temperatures, but it fails at the extremely high 
temperatures encountered in hypersonic flight. The viscosity begins deviating 
from the Sutherland value due to the onset of dissociation at approximately 
3000 K. The viscosity near the wall is also affected by the extremely low 
pressures encountered at high altitude [11] . This effect is due to velocity slip 
at the wall (see next section) . The Sutherland formula, equation (47) , is 
recommended for temperatures less than 3000 K. 




su 


1. 458 x 10~ 5 T 3/2 
110.4 + T 


[gm/cm sec] 


(47) 


This formula is_also used in the computer program (Appendix A) to 
calculate using T^ computed from the adiabatic relation for temperature, 
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although it is realized that neither g nor has valid physical meaning at 

extremely high stagnation enthalpy. This use is justified because Ji is used 
only for nondimensionalizing jl and for computing Re . 

you 

The viscosity of dissociated and ionized air is approximated in Reference 
6 using a simple summation formula for a mixture of hard spherical molecules: 


JL_ 




su 




(48) 


where 

ju gu = viscosity at same temperature from Sutherland formula 

W = equivalent molecular weight of undissociated air 
K 

= mean free path of species i 
= mean free path of undissociated air molecules. 


The ratio of mean free paths in equation (48) is given by 
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where 


f 


S.^ = collision cross section for particle i with particle j 

S = collision cross section for undissociated air molecules. 
R 


(49) 


The collision cross sections are tabulated as a function of temperature in 
Reference 6. 
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The viscosity given by equation (48) is a function of the composition of 
the gas mixture and also of the temperature via the dependence of S R and S„ 

on temperature. Since in the present analysis the air is not in chemical equi- 
librium, the viscosity obtained from the computer program at a given tempera- 
ture differs greatly from that shown in Reference 6 for equilibrium conditions. 
As the mass fractions of the components of air approach their undissociated 
values, the numerical value of the viscosity ratio in equation (48) approaches 
1.0, i.e., n approaches the Sutherland value as expected. 

3. Boundary Conditions . To solve the governing equations given in the 
previous section, freestream and wall boundary conditions are required. 

Freestream (rj = 1): 

The air at the freestream boundary is in its undisturbed state. 
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Wall (tj = 0); 


In high -altitude, low-Reynolds number flight, the continuum model of the 
gas breaks down in regions of large gradients of the physical properties near the 
wall. Hence, the Navier-Stokes description is invalid for the gas layer near the 
wall (Knudsen layer) with thickness on the order of the mean free path [11] . 

Also, the familiar continuum zero-velocity and zero -temperature -jump wall 
boundary conditions are not applicable. Although the Navier-Stokes equations 
are invalid near the wall, they can still be used, down to quite low Reynolds 
numbers, to describe the outer flow field if the proper boundary conditions are 
used at the outer edge of the Knudsen layer. These boundary conditions, known 
as slip conditions, are the mean velocity, temperature, and species mass 
fractions. To compute these slip conditions, a kinetic theory approach must be 
used for the Knudsen layer. The boundary conditions for the Knudsen layer are 
the mean slip conditions at the outer edge and the kinetic gas-surface conditions 
at the wall. 

The Boltzmann equation is the governing equation for the kinetic theory 
description of a flow field. However, due to the difficulty in solving the Boltz- 
mann equation for the Knudsen layer, we resorted to an approximate kinetic 
theory slip model. Reference 4 gives the details of the derivation of the slip 
conditions for a multicomponent reacting gas. By matching the species, momenta, 
and energy fluxes at the outer edge of the Knudsen layer to the difference between 
the incident and reflected fluxes at the wall, the jump in the desired properties 
across the Knudsen layer is obtained. The fluxes are calculated by taking 
moments of the velocity distribution function which is approximated by using a 
Chapman-Enskog expansion for a multicomponent mixture. The species flux is 
greatly affected by the catalytic nature of the wall. The wall is assumed to be 
catalytic with respect to recombination of dissociated molecules. Equations 
are obtained for a partially catalytic wall. The extremes of noncatalytic and 
fully catalytic walls are easily obtained from the equations for a partially 
catalytic wall. The resulting nondime ns ional equations for slip velocity, 
temperature, and species are as follows: 
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( 50 ) 



s 

The superscript, s, on Y. denotes the value of Y. 
Knudsen layer. 1 


at the outer edge of the 
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The similarity equations (18) through (24) and the coordinate trans- 
formation equation (25) are used in equations (50) through (52) to yield the 
following slip equations: 



( 54 ) 


and 
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The quantity y in equation (55) is the recombination coefficient for the 
i species; it is a measure of the catalyticity of the wall. The value of varies 
from 0 (noncatalytic wall) to 1 (fully catalytic wall). 

The computer program used in this report has the fully catalytic wall and 
noncatalytic wall options, but it does not, in the present form, have a partially 

g 

catalytic wall option [for a given catalyticity y., equation (55) would give Y. ]. 

In the fully catalytic wall option, the wall is assumed to be catalytic only with 
respect to recombination of neutral atomic species; it is assumed to be non- 
catalytic with respect to recombination of the charged particles NO and e“. 

The wall boundary conditions on Y. for these special cases are given as 
follows: 

Noncatalytic Wall: 

For a noncatalytic surface (y = 0), equation (55) reduces to 


Z 


(y;.) = 

1] S 


(56) 


A sufficient condition for equation (56) to be satisfied is that 


(Yjj) s =0 j = 1 NS . (57) 


Equation (57) is also the most physically plausible means by which equation (56) 
can be satisfied; therefore, equation (57) is taken as the noncatalytic boundary 
condition for species mass fractions. 
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Fully Catalytic Wall: 


For the fully catalytic wall, the surface is assumed to be fully catalytic 
with respect to recombination of neutral atomic species, but it is assumed to be 
noncatalytic with respect to recombination of the charged particles NO + and e“. 
Therefore, the effect of the wall is to drive the gas towards its freestream 
composition, except for the charged particles NO + and e~. The following 
boundary conditions are then obtained: 


( Y »«). ■ 

0.767 

(58) 

M. ■ 

0.233 

(59) 

^ Y NO^s 

(Vs - Vs ■ 0 

(60) 

( Y NO+)s 

V>s - 0 • 

(61) 


Equations (58), (59), and (61) produce the impossible result that 
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however, because (Y. T _. ) and (Y _) 

. . , . v NO +/ s v e i 

tolerable. 


are usually very small, this error is 


B. Method of Solution 

Equations (26) through (32) and (33) through (39) together with equations 
(46) and (48) constitute two sets of nonlinear, coupled ordinary differential equa- 
tions with boundary conditions previously given. The first order equations are 
solved by direct numerical quadrature, and the second order equations are 
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integrated by a finite difference method known as the successive accelerated 
replacement method which is an iterative scheme that starts from a guessed 
solution. The salient feature of this method, proposed by Lieberstein [13], is 
that the successive corrections applied to the flow profiles in each iteration are 
controlled by acceleration factors which are used to increase the rate of con- 
vergence of the computed quantities. Thus, this method can be successfully 
applied even if the initial, guessed profiles do not approximate the converged 
solutions very well. For the present application, this statement holds true in 
the midrange of the flow regime for which our analysis is applicable. However, 
divergence problems are encountered at the continuum end of the regime 
(Appendix A) . 


C. Computer Program 

Details of the computer program are found in Reference 10, while 
Appendix A is a current users manual. 


III. RESULTS 


Results from the computer program, with the modifications introduced 
in this report, are first compared with the ’’old" program. Then, the program 
output is compared with some available experimental data. Finally, some 
general results from the program and input data for running the program are 
presented. 


A. Effect of Computer Program Modifications 

The program described as the old program includes all the analysis in 
this report except the constitutive equations (46) and (48) for enthalpy and 
viscosity, respectively. The old enthalpy equation assumed a fully excited 
vibrational state for diatomic molecules. This assumption is always violated 
near the freestream and near the wall because temperatures in these regions 
are less than ® The old viscosity equation (Sutherland) is restricted to non- 

dissociated air or to temperatures less than approximately 3000 K at equi- 
librium conditions. This temperature is greatly exceeded in the central part 
of the shock layer at reentry speeds. 
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B. Comparison of Data Between Old Program and 
Modified Program 

Computed values of shock layer thickness, n ; wall heat transfer 
coefficient at the stagnation point, C„; and maximum ratio of viscosity to 

XI 

Sutherland viscosity on the stagnation streamline, p/p , are given in Table 1 

for a 30.5-cm radius sphere using (l) the ”old M program, (2) the program with 
enthalpy modification only, and (3) the program with enthalpy and viscosity 
modifications. The conditions for which the runs were made are listed in 
Table 1. The thermodynamic properties associated with altitude are obtained 
throughout this report from Reference 14. 

TABLE 1. DATA COMPARISON 

alt = 86 km 


r, = 30. 5 cm 
b 

V = 793 000 cm/sec 

00 


Program Identification 

Shock 

Layer 

Thickness 

n 

00 

Wall Heat 
Transfer 
Coefficient 

C H 

(f~) 

\ sv /max 

1. Old Program 

0. 1361 

0.147 

1.0 

2. Modified Program 

(a) New enthalpy computation, 
equation (46) 

(b) Sutherland viscosity 

0.1409 

0.154 

1.0 

3. Modified Program 

(a) New enthalpy computation, 
equation (46) 

(b) New viscosity, 
equation (48) 

0.1411 

0.172 

1.37 
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For these runs, conditions were chosen at which the modifications in 
viscosity and enthalpy computation would have their greatest effect. These 
conditions are the ones which yield a high degree of dissociation and ionization, 
equations (46) and (48), i.e. high speed and high freestream density. There- 
fore, an altitude was chosen (86 km) where freestream density is near the 
maximum for which the program will run (see General Results section) . The 
table shows that the enthalpy and viscosity modifications had a significant effect 
on heat transfer rate and shock layer thickness. The heat transfer rate was 
increased 17 percent due to the combined effect of both modifications. 


C. Comparison with Experimental Data 

Very little experimental data exist for the low density, hypersonic, high 
stagnation enthalpy flow regime for which this computer program was designed. 
Therefore, we do not have the detailed comparison with experimental data which 
is desired to validate the present analysis. However, a limited amount of heat 
transfer and pressure data for spheres are available from an arc-jet facility at 
moderate stagnation enthalpy [15] , and a larger body of stagnation heat transfer 
data exists for relatively low enthalpy, and hypersonic flow [16] . 

Figure 2 compares computed stagnation point heat transfer coefficients, 

C , for a sphere with experimental data from References 15 and 16. The data 
H 

in this figure represent several levels of stagnation enthalpy and, hence, several 
levels of real gas effects, so good agreement among all the data is not expected. 
The purpose here is to compare the computed results directly with each set of 
experimental data and to observe the trend in C with changing stagnation 

enthalpy. A best fit curve through the relatively low enthalpy data of Vidal and 
Wittliff [ 16 ] is shown by the solid line , and the square symbols give the computed 
values for the same flow conditions. Although the computer program was 
designed for real gas, chemically reacting flows, it can be run at the low 
enthalpy conditions of the Vidal and Wittliff data, in which case the computed 
dissociation levels are very low. The data are plotted as a function of the 
similarity parameter, K 2 , used by Vidal and Wittliff, and defined as follows: 
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Figure 2. Comparison of present theory with experimental data for stagnation 

point heating on a sphere. 
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(63) 


and 


T* - | <V + V , (64) 

where 

H = viscosity computed by Sutherland’s formula using T* 

= viscosity computed by Sutherland’s formula using . 


The computed data do not agree with the mean of the Vidal and Wittliff 
data as well as we would like, but the computed data points fall within the 
scatter of the experimental points [16], The computed values tend to over- 
estimate C at rarefied conditions and underestimate it for denser flow. It is 
H 

important to note that at the higher enthalpy cases, the wall catalysis can affect 

the value of C by a factor of 3 or more [4], Hence, any comparison of experi- 
H 

mental and calculated heating must be for similar wall catalysis. 

The arc-jet data of Scott [15] are shown by the filled diamond symbols, 
and computed data for the same conditions are shown by the empty diamond 
symbols. The flow conditions for these tests are given in Table 2. 

The free stream boundary conditions in the computer program were 
modified to match the dissociated arc-jet freestream. Also, the spheres in the 
arc-jet tests were coated with teflon to produce a noncatalytic wall and, there- 
fore, the noncatalytic wall option in the computer program was used for 
comparison with these data points. The noncatalytic wall option was also used 
for the remainder of the computer data in this figure and the entire report 
unless otherwise noted. 

The comparison of the computed data with Scott’s experimental data is 
considered acceptable. 
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TABLE 2. ARC-JET FLOW FIELD PROPERTIES 



Case 1 

Mach number, 

10 

Reynolds number, Re ro 

380 

Freestream speed, V (cm/sec) 

423 000 

Freestream temperature, T^ (K) 

270 

Freestream density, (gm/cm 3 ) 

3. 28 x 10- 8 

Wall temperature, (K) 

450 

Radius of sphere, r^ (cm) 

5.08 

Species mass fractions in freestream: 


n 2 

0.4138 

o 2 

0.6519 x 10- 5 

N 

0.3514 

O 

0.2348 

NO 

0.5380 x 10 -5 

NO + 

0.2221 x 10 “ 5 

e” 

0.3850 x 10“ 10 


Case 2 

8 

553 

412 000 
414 

6.41 x 10 -8 

450 

5.08 

0.6124 

0. 1084 x 10" 4 
0. 1522 
0.2347 

0. 3761 x 10 -4 
0.6849 x 10 -5 
0.1187 x 10 -9 


Some high-enthalpy computed data, for which the dissociation and ioniza- 
tion levels are very high, are shown in Figure 2 by the circular symbols. These 
data points fall considerably below the lower enthalpy data at the continuum end 
of the flow regime (at large values of K 2 ) where there is a consistent trend of 

decreasing C TT with increasing enthalpy for a given K 2 . The parameter K 2 does 
H 

not account for real gas effects and, therefore, should not be expected to corre- 
late data with widely different enthalpy levels. 
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Acceptable comparisons have also been obtained with other theoretical 
predictions by Kumar and Jain [3] who used an analysis identical to the present 
one except for slip, enthalpy, and viscosity computations. 


D. General Results 

Figures 3 and 4 illustrate typical convergence behavior of the solutions 
with increasing number of iterations. Figure 3 shows the successive computed 
nondimensional temperature profiles, including the initial guess, along the stag- 
nation streamline of a sphere with noncatalytic wall. The temperature profiles 
did not quite converge by the final iteration (2000) , so it was necessary to make 
a rerun using the output of run No. 1 as input data for run No. 2. The figure 
illustrates that extremely accurate initial guesses are not generally required to 
achieve convergence; however, the more accurate the initial guess, the faster 
the solutions will converge. At very large Reynolds numbers, where gradients 
in flow properties are large, more accurate initial guesses are required to 
avoid divergence. Figure 4 shows the convergence behavior of wall heat trans- 
fer coefficient, C TT , for the same run. The converged value of C„, obtained in 
H ri 

run No. 2, is shown at the right margin. 

Figure 5a gives the computed velocity components, and the pressure, 
temperature, and density profiles along the stagnation streamline of a sphere 
with noncatalytic wall at relatively high freestream Reynolds number (1315). 
Figure 5b gives the mass fraction profiles for the same run. This rim is 
approaching the maximum Reynolds number for which the program will converge 
easily without extremely accurate input data or modifications to the program such 
as spacing the computation increments closer together. The meaning of the 
quantities in Figure 5a are given in equations (18) through (24) ; e.g. P t is the 
pressure on the stagnation streamline, and P 2 gives the correction in pressure 
for small angles away from the stagnation streamline. A thick ’’shock wave” is 
indicated by the steep gradients in the region 1.06 < r < 1. 12. 

The same quantities presented in Figure 5 are plotted in Figure 6 for a 
smaller Reynolds number (185) . Figure 6a reveals that all evidence of a shock 
wave has disappeared at this low Reynolds number. Comparison of Figures 5b 
and 6b shows that the degree of dissociation and ionization decreases drastically 
in lower density flow. 

Figure 7 gives the flow profiles along the stagnation streamline of a 
sphere with fully catalytic wall at a medium Reynolds number (458) . The 
Sutherland viscosity formula was used because this rim was made before the 
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Figure 4. Typical convergence of stagnation point heat transfer coefficient for sphere with 

increasing number of iterations. 




a. Velocity and thermodynamic property profiles. 

Figure 5. Flow profiles for sphere with noncatalytic wall at large Reynolds number (Re ro = 1315) 
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b. Mass fraction profiles. 
Figure 5. (Concluded) 
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a. Velocity and thermodynamic property profiles. 


Figure 6. Flow profiles for sphere with noncatalytic wall at small Reynolds number (Re co = 185). 
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viscosity modification to the program. The density profile gives a hint of a 
shock wave similar to that shown in Figure 5a. Comparison of Figures 5b, 6b, 
and 7b illustrate the difference in species mass fraction profiles between a 
noncatalytic wall and a fully catalytic wall. The species mass fractions approach 
their freestream values at the catalytic wall except for the ions NO + and e“. The 
wall is assumed noncatalytic with respect to ion recombination in this program. 

Figures 5 through 7 provide useful starting data for running the program 
at conditions near the ones in these figures. 

Figure 8 compares temperature, density, and velocity profiles for a 
sphere and cylinder at the same flow conditions as in Figure 5. The profiles 
are remarkably similar. Therefore, the profiles for a sphere, Figures 5, 6, 
and 7, are adequate to use as starting data for the cylindrical option of the 
program. 

The change in the stagnation line temperature profile with freestream 
density, or Re^ , is illustrated in Figure 9. The peak value in nondimensional 

temperature decreases and the shock layer thickness increases with increasing 
altitude or decreasing density. 

The freestream density and speed greatly affect the degree of dissociation 
of the air molecules and, hence, the viscosity of the gas mixture. Figure 10 gives 
the peak values from the viscosity and mass fraction profiles for N, O, and NO 
as a function of altitude for a fixed large value of freestream speed. As before, 
the freestream thermodynamic properties associated with altitude are obtained 
from Reference 14. Although the stagnation enthalpy is large, the figure shows 
that dissociation becomes negligible and viscosity approaches the Sutherland 
value at altitudes greater than approximately 100 km. This is due to the 
decreased reaction rates at the lower temperatures which occur in the shock 
layer at higher altitudes (Fig, 9). Figure 11 illustrates that dissociation and 
the viscosity ratio increase with increasing freestream speed at the fixed altitude 
of 96 km. 

A parametric study was made to determine the effect of altitude (or 

density) , freestream speed, and wall temperature on the stagnation point heat 

transfer coefficient, C TT . Figure 12 gives C TT for a sphere and cylinder as a 

H xl 

function of altitude with fixed freestream speed and wall temperature. Hansen’s 
viscosity, equation (48), was used for the sphere, but Sutherland’s viscosity 
was used for the cylinder because those runs were made before the program 
was modified to compute Hansen's viscosity. The figure shows a steady increase 
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Figure 9. Temperature profiles at various Reynolds numbers for sphere 
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Figure 10. Effect of variation in altitude on dissociation and viscosity 

in shock layer of sphere. 
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Figure 11. Effect of freestream speed on dissociation and viscosity 

in shock layer of sphere. 
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Figure 12. Stagnation point heat transfer coefficients for spheres and 
cylinders as a function of altitude. 


in C with increasing altitude. The theoretical free molecular value (C = 1.0) 

is exceeded for the sphere at approximately 106 km and for the cylinder at 
approximately 109 km. The curves continue to climb with increasing slope at 
higher altitudes. This unrealistic behavior for a fixed wall temperature might 
be attributed to a breakdown in the continuum approach (Navier-Stokes equations 
with slip boundary conditions) at extremely high altitudes and low densities. 

In particular, the assumption of a thin Knudsen layer, which is implicit in the 
present analysis, becomes invalid at very low densities. One must also consider 
the practicality of a specified T^ in rarefied flow. This breakdown is gradual, 

and one cannot pinpoint a sharp boundary beyond which the program gives invalid 
results. However, the results should be treated with increasing skepticism, for 
this body size and flow conditions, at altitudes greater than approximately 

104 km, or He < 40. 

00 

An interesting feature of Figure 12 is the crossover of the curves for 

the sphere and cylinder. This crossover should not be attributed to the differences 

in viscosity computation because the effects of Hansen’s viscosity (used for the 

sphere but not for the cylinder) is to increase C at low altitudes but not to 

fl 

affect C TT at high altitudes where dissociation is negligible. Therefore, if 

Hansen’s viscosity were also used for the cylinder, the crossover of the two 
curves should be expected at a higher altitude. The explanation for this cross- 
over is a result of the stronger merged shock layer effect on the sphere. 

The variation of C TT with freestream speed for a sphere at fixed altitude 
ii 

(96 km) and wall temperature (600 K) is shown in Figure 13. The effect of 
increasing freestream speed, or stagnation enthalpy, is to decrease C . 

Figure 14 gives the variation in C with wall temperature for a sphere 

H 

at fixed altitude (96 km) and speed (793 000 cm/sec). An increase in wall 

temperature produces an increase in C . 

rl 

The data in Figures 12 through 14 are replotted in Figure 15 as a function 

of K 2 ; the similarity parameter is given by equation (62). This parameter 

includes all the independent variables in Figures 12, 13, and 14, but the effect 

of freestream density is obviously dominant in determining the value of K 2 and 

C . The data were also plotted as a function of other similarity parameters 
H 

(not shown) , but K 2 correlated the data as well as, or better than, any of the 
other parameters. This is not the case if wall catalysis is varied. Therefore, 

K 2 was selected as the similarity parameter to use in presenting the remainder 
of the data. 
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Figure 15. Stagnation point heat transfer coefficient for sphere 

as function of K 2 . 
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Figure 16 gives the nondimensional shock layer thickness, n m , for 

spheres and cylinders as a function of K 2 . The shock layer is thin at the con- 
tinuum end (large K 2 ) of the flow regime and becomes very thick at highly 
rarefied conditions (small K 2 ). 

Figure 17 presents the computed slip speed for a sphere and cylinder as 
a function of K 2 . The slip speed is computed at the outer edge of the Knudsen 
layer, but it is applied at the wall in the outer flow solution, i.e. , a thin Knudsen 
layer is assumed. This is a source of error at high altitudes where the Knudsen 
layer thickness becomes significant. The slip speed is greatly affected by varia- 
tion in freestream speed and wall temperature. However, these effects are not 
correlated well by the parameter K 2 as shown for the sphere by the partially 
filled symbols. To estimate the slip speed for given conditions, one can 
extrapolate from the solid curves using the partially filled symbols as a guide in 
correcting for freestream speed and wall temperature. Once again it appears 

that the proper T must be used for rarefied flow calculations. The decrease 

w 

in slip speed with increasing altitude (decreasing K 2 ) at high altitude (small K 2 ) 

seems unrealistic and has been criticized. However, decreasing slip speed 

decreases the C computed by this program. Since the computed C TT is too 
xi H 

large at high altitudes for a fixed T^ (Figs. 12 and 15) , the decrease in slip 

speed at high altitudes tends to compensate for the overprediction of C and is, 
from a practical standpoint, a favorable phenomenon. 

The temperature at the edge of the Knudsen layer is given as a function of 
K 2 in Figure 18. This temperature also depends strongly on freestream speed 
and wall temperature. 

The information in Figures 15 through 18 is useful as input data in 
running the computer program (Appendix A) . 

An example of the practical use of this computer program is given in 

Figure 19. Stagnation point heating rates for the spherical nose of the Space 

Shuttle External Tank are given as a function of trajectory time. The solid 

curve is the heating rate predicted by the MSFC Thermal Environment Branch 

(ED33) and is based on theoretical continuum methods at low altitudes and 

experimental data at high altitudes. The filled diamond symbols were computed 

using the present program. The trajectory time interval from approximately 

275 to 475 sec corresponds to an altitude interval which is too high for this 

program, i.e. , the predicted values of C„ were unrealistically large. For a 

±1 
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6, HEATING RATE 





body of this size (r^ = 30. 5 cm) , the altitude interval for which this program is 

applicable is approximately 86 to 105 km. For this altitude interval, the pre- 
dicted values from this program agree quite well with the predictions from 
experimental data for a known T . 

w 
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APPENDIX A. COMPUTER PROGRAM 



A P PEN D I X A. COM PUTER PROG RAM 


General Information 


The computer program is written in Fortran IV language , and it is run 
on the Univac 1108 computer at MSFC. As currently written, the program is 
limited to a maximum of 2000 iterations. If convergence is not achieved in 
2000 iterations, the program is rerun using the output from the first run as 
starting data for the second run. 

The program has three options: 

(1) Body: sphere or two-dimensional circular cylinder 

(2) Wall catalyticity: noncatalytic wall or fully catalytic wall 

(3) Viscosity: Sutherland or Hansen's high temperature model. 

The choices made in options (1) and (2) make very little difference in 
the computer time used. However, the Hansen viscosity option takes significantly 
more computer time than the Sutherland option; e.g. , if the full 2000 iterations 
are used, the Sutherland option takes approxi mately 8 min and the Hansen option 
takes about 10 min. 


Program Input 


The quantities needed for input to the program are defined and their 
functions described as follows: 


Computer 

Program 

Variable 

Symbol 
Used in 
Report 

Dimensions 

Description 

Suggested 
Source or 
Range of Values 

EFFR1 

n 

CO 

Dimensionless 

Initial guess for nondime ns ion al 
shock layer thickness 

Figure 16 

TOL 


Dimensionless 

Computation is stopped when 

C u converges within a pre- 
ri 

scribed tolerance. If 
IC -C I < TOL , 

1 h n Vi 1 

stop program and print results 

0.00001 - 0.0001 

EPSI 


Dimensionless 

Maximum change allowed in all 
computed dimensional quanti- 
ties from one iteration to the 
next 

0.001 - 0.005 
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Computer 

Program 

Variable 

Symbol 
Used in 
Report 

Dimensions 

Description 

Suggested 
Source or 
Range of Values 

APT 


km 

Altitude corresponding to the 
free stream thermodynamic 
properties (used for identifica- 
tion only) 


RADB 

? b 

cm 

Radius of body 


UFS 


cm/s 

Free stream speed 


TFS 


K 

Free stream temperature 


TWK 

T 

W 

K 

Temperature of body surface 


ROFS 

Pa, 

gm/cm 3 

Free stream density 


U(l) 

u = u /sine 
Iw W 

Dimensionless 

Initial guess for (slip speed) /sin# 

Figure 17 

T(l) 

T i 

lw 

Dimensionless 

Initial guess for nondime ns ional 
temperature of gas at outer 
edge of Knudsen layer 

Figure 18 

U(I) 

U 1 

Dimensionless 

Initial guess for u^ (see equa- 
tion 18) at location I (I = 

2, 3, 50) 

Figures 5-7 

V(I) 

V 1 

Dimensionless 

Initial guess for (see equa- 

tion 19) at location 1(1 = 

2, 3, 50) 

Figures 5-7 

T(l) 

T 1 

Dimensionless 

Initial guess for T^ (see equa- 
tion 20) at location I (I = 

2, 3, •••» 50) 

Figures 5-7 

RO(I) 

P 1 

Dimensionless 

Initial guess for (see equa- 

tion 21) at location I (I = 

2, 3, 50) 

Figures 5-7 

P2(I) 

P 2 

Dimensionless 

Initial guess for P 2 (see equa- 
tion 22) at location I (I = 

2, 3, . . . , 50) 

Figures 5-7 

C02(I) 

O* 

to 

Dimensionless 

Initial guess for mass fraction 
of 0 2 at location I (I = 

2, 3, 50) 

Figures 5-7 




Computer 

Program 

Variable 

Symbol 
Used in 

Dimensions 

Description 

Suggested 
Source or 
R ange of V alue s 

CNOI(I) 

y no + 

Dimensionless 

Initial guess for mass fraction 
of NO + at location 1(1 = 

2 f 3, < < • j 50) 

Figures 5-7 

CNO(I) 

y no 

Dimensionless 

Initial guess for mass fraction 
of NO at location I (I = 

2 f 3, • ■ • f 50) 

Figures 5-7 

CN(I) 

y n 

Dimensionless 

Initial guess for mass fraction 
of N at location I (I = 

2 9 3, • • • f 50) 

Figures 5-7 

CO(I) 

Y o 

Dimensionless 

Initial guess for mass fraction 
of O at location 1(1 = 

2 f 3 f « • • y 50) 

Figures 5-7 


The quantities U(I) through CO(I) must be input for each of the 49 
equally-spaced points in the shock layer along the stagnation streamline starting 
with 1=2 near the wall and ending with I = 50 near the freestream (see Fig. A-l) . 
The wall (1=1) and freestream (I = 51) boundary conditions are given by the 
quantities UFS through T(l) . 


The format for the input data is as follows: 


Card 

Format 

Variable 

Name(s) 

Location and Description 

1 

12A6 

AB1, . .., 
AB/2 

Columns 1-72 contain a comment for 
project identification 

2 

3F10.0 

XBOD 

Columns 1-10 contain the body option 
(sphere = 1.0, cylinder = 2.0) 



XCAT 

Columns 11-20 contain the wall 
catalyticity option (noncatalytic wall = 
1.0, fully catalytic wall) 



XVIS 

Columns 21-30 contain the viscosity 
option (Sutherland viscosity = 1.0, 
Hansen viscosity = 2.0) 

3 

2F10.0 

CASE 

Columns 1-10 contain the run number. 
This identification is useful to distinguish 
between runs whenever a rerun is 
necessary 


62 



Card 

Format 

Variable 

Name(s) 

Location and Description 



EFFR1 

Columns 11-20 

4 

2F10.0 

TOL 

Columns 1-10 



EPSI 

Columns 11-20 

5 

5F10.0 

ALT 

Columns 1-10 



RADB 

Columns 11-20 



UFS 

Columns 21-30 



TFS 

Columns 31-40 



TWK 

Columns 41-50 

6 

E13.8 

ROFS 

Columns 1-13 

7 

2F10.0 

U(l) 

Columns 1-10 



T( 1) 

Columns 11-20 

8-56 

5E16.9 

U(I) 

Columns 1-16 



V(I) 

Columns 17-32 



T(I) 

Columns 33-48 



RO(I) 

Columns 49-64 



P2(I) 

Columns 65-80 

57-105 

5E16.9 

C02(I) 

Columns 1-16 



CNOI(i) 

Columns 17-32 



CNO(I) 

Columns 33-48 



CN(I) 

Columns 49-64 



CO(I) 

Columns 65-80 
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PROGRAM OUTPUT 


The program printout is described in this section. Figure 5 of the text 
was obtained from this printout. Explanation of the printout follows. 


Title Page 

Options used — The options which were selected for body, wall cata- 
lyticity, and viscosity are printed out here. 

Input data — The input data are printed out for record. 

Computed freestream data — Some freestream quantities which are used 
for nondimensionalizing and for computation of various similarity parameters 
are printed here. The equations used are as follows: 


VSOUND = n/ y R T ffl [cm/sec] 


(A-l) 


where 


y = 1.4 

R = 2. 8708 x 10 6 cm 2 /sec 2 K . 


M = AMACH 

CO 


UFS 

VSOUND 


(A-2) 


T 


Qco 


STAGFS = T 

CO 



[ °K] 


1.458 x 10" 5 T^ 

= RTFS = — [gm/cm-sec] 

T + H0.4 


1.458 x 10" 5 f 

- 0°3 

p n = RTSTAG = [gm/cm*sec] . 

V + n0 - 4 


(A-3) 
(A -4) 


(A-5) 
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It is recognized that the stagnation temperature, equation (A-3), has 
little physical meaning at the high stagnation enthalpies for which this computer 
program is intended. This is because equation (A-3) is based on the assumption 
of constant specific heat for temperatures up to ; the assumption of constant 

specific heat is violated at moderate temperatures. However, is used only 

for nondimensionalizing and for computing . It does not influence the compu- 
tation of the flow field in any way. Likewise, p is a fictitious viscosity 

because T. is fictitious and also because the Sutherland formula is not applic- 
0°° 

able at high temperatures. 

Similarity Parameters — Some similarity parameters are printed out 
here. Several of these have been used, with some success, in correlating 
hypersonic or low density data. The similarity parameters are defined as 
follows: 


^co v co"b , 

Re = RE YF = = , (A-6) 

°° u 

P CO 


p V r, 

H CO CO b 

Re = REYN = z 

000 


(A-7) 


p V r, 

; co CO b 

Re =s REWALL = = 

w a 

'w 


(A-8) 


where ^ is computed by the Sutherland formula using T^ , 


K = XKNFS = — , 

n r. 

co b 


(A-9) 
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(A -10) 



Data Pages 


A set of computed data is printed out after every 200 iterations. Each 
set of data consists of three pages. 

1. First Page . The computed shock layer thickness (EFFE) and wall 

heat transfer coefficient (C ) are printed at the top of the page after every 50 

H 

iterations. The control factors which are used to update the computed quantities 
to obtain input variables for the next iteration are also printed here. 

The columns on this page give nondimensional values of ( 1) radial posi- 
tion from the center of the body to the equally-spaced points at which the inde- 
pendent variables are computed, (2) the velocity component Uj [see equation 
(18)], (3) the velocity component Vj [see equation (19) ] , (4) specific enthalpy, 
(5) temperature, (6) viscosity ratio, local ratio of Hansen viscosity to Suther- 
land viscosity, (7) pressure Pj [see equation (22)], (8) pressure correction 
Pj [see equation (22)], and (9) density. 

Note that VFS is the same as UFS elsewhere in the program. It is the 
free stream speed . 

2. Second Page . This page gives the reaction rates for each of the 
seven species. 

3. Third Page . This page gives the mass fractions of the chemical 
species which compose the air. Note that mass fraction was denoted by ”Y M 
in the text rather than by M C.” The electron number density is also given. 

When C converges within the prescribed tolerance, or when the 
H 

iterations reach 2000, the iterations are stopped and the final data set is 
printed out. At this point a set of 98 cards is punched to record the velocity 
and thermodynamic data given on page 1 and the mass fraction data on page 3 
of the data printout. This set of cards can then be used as starting data for 
a new run. 

To determine whether or not the program has converged sufficiently and 
whether the results are reasonable, the following checks are suggested. 

(1) Check EFFE for convergence 

(2) Check C for convergence 

H 
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(3) Check the slip conditions for convergence, i.e. compare U1 and T 
at RAD/RADB =1.0 for the last few sets of data 

(4) Scan the velocity component and thermodynamic property columns 
to determine if the profiles are smooth and the values are reasonable 

(5) Scan the mass fraction columns to check for reasonable behavior. 
Particularly check the N 2 column near the wall to determine if these mass 
fractions exceed the free stream value. 


SOME DIFFICULTIES ENCOUNTERED IN RUNNING PROGRAM 


The program generally runs without much difficulty if the starting data 
are within reasonable bounds; however, problems do occur, and some of the 
most common problems are discussed as follows. 

1. Failure of EFFR to Converge . The shock layer thickness, EFFR, 
usually is the first quantity to converge; however, it sometimes diverges. The 
most common cause of this is using starting data from a run with a much 
different stagnation enthalpy, or T^ , from the case being run. This causes 

a discontinuity in the nondimensional temperature profile between the points 
I = 50 and I = 51 because the point I = 50 comes from the starting data from a 
previous run and point I = 51 is the free stream boundary condition. The 
decision on whether to decrease, increase, or not change EFFR is made on 
the basis of the temperature gradient near the freestream (see Appendix B) . 
Therefore, if a sizable discontinuity in temperature exists at this location, 
EFFR will diverge. 

One solution to this problem is to modify the last few input cards to make 
a smooth transition in the temperature profile to the freestream temperature. 
Probably a_better solution would be to change the nondimensionalizing tempera- 
ture from T~ to T . This would guarantee a smooth transition, in the 

0°0 03 

starting data, to the freestream temperature. To prevent divergence of EFFR, 
the program has been modified to restrict EFFR within the range 
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\ EFFR1 < EFFR < f EFFR1 

/ A 


If EFFR1 is not guessed closely enough, this restriction will necessitate a 
rerun using a different EFFR1. 

2. Unrealistic Chemical Composition . If the species mass fractions in 
the starting data are greatly different from the true values for the case being 
run, unrealistic computed mass fractions can result. This can lead to unrealis- 
tic computed values for all other quantities. For example, this might occur for 
a rather low speed, high density case in which dissociation should be negligible 
(Fig. 11) . However, since starting data from a low speed, high density run 
were not available, data were used from a high speed, high density run in which 
the degree of dissociation was high. The program could not handle the grossly 
incorrect starting data, and the computed results became more and more 
unrealistic from one iteration to the next. The unrealistic data included mass 

fractions for N 2 much larger than the freestream value, negative C (heat 

H 

transferred from the body to the gas) and a monotonically increasing EFFR 
beyond reasonable bounds (Fig. 16). Note that was called "STANC" in the 
older printout. 


SOME SUGGESTIONS FOR RUNNING PROGRAM 


To minimize the previously discussed problems and to facilitate operation 
of the program, the following suggestions on input data are made. 

1. EFFR1, U(l), T(l), and TOL. Figures 16, 17, and 18 can be used 
for obtaining initial guesses for EFFR1, U(l), andT(l), respectively. The 
value for TOL probably should be set smaller than the desired tolerance on C 

H 

to assure that the other quantities converge also. A value of 0.00001 for TOL 
was used in most of the runs made for this report. 

2. Velocity, Thermodynamic Property, and Mass Fraction Profiles . 

For best results, the starting data profiles for the quantities U(I) through 
CO(I) should be obtained from a rim which matches the required freestream 
conditions of the new run as closely as possible. Close matching is especially 
important at low altitude, high density conditions for which gradients in the 
flow properties are large. The question arises as to what constitutes a 
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sufficiently close ’’match” and which matching, or similarity, parameter should 
be used. If every free stream and body condition is matched fairly closely, it is 
unnecessary to use a similarity parameter. From experience, it is found that, 
for a given body, changes in freestream density, p ro , and speed, V , in that 

order, have the greatest effect on the outcome of the run. Therefore, for 
bodies of the same size, one should try to match and V m . A good rule of 

thumb is that p^ for the new rim should not differ from p m for the starting 

data run by more than the equivalent of 3 km altitude in the standard atmosphere. 
Near the low altitude extreme of the program’s capability, this difference should 
not be more than approximately 1 km. A difference in of approximately 

2000 m/sec can usually be tolerated. For bodies of greatly different size, one 
should probably use a similarity parameter, such as K 2 , to match starting data 
to the new run. 

Figures 5, 6, and 7 can be used for starting data. If the desired free- 
stream properties for the new run are farther removed from those of Figures 5, 
6, and 7 than the previously suggested increments, it is advisable to reach the 
desired conditions in two or more runs. If severe difficulties like those pre- 
viously discussed are encountered, it is probably best to go back to a ’’good, ” 
or converged, set of starting data and use smaller increments in freestream 
properties rather than to use the output of the ”bad” run as starting data for a 
rerun. If the program is used frequently, a library of starting data card sets 
can be built up to cover the range of possible freestream conditions for which 
the program is applicable. 

3. Hansen-Sutherland Viscosity Option . If the computed ratio of 
Hansen’s viscosity to Sutherland’s viscosity is approximately 1.0 for a given 
set of freestream conditions, one can use the Sutherland viscosity option for 
nearby freestream conditions to reduce computer time. 
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APPENDIX B. LISTING OF PROGRAM 


I 



<1 


MAIN 


DATE 120877 


0C1CC 

1* 

aoioo 

2* 

CC100 

3* 

00130 

4 * 

00100 

5* 

00100 

5* 

CP10C 

7* 

00101 

8* 

C0101 

9* 

30101 

10* 

00103 

11* 

CC103 

12* 

0C1D4 

13* 

001C5 

14 * 

CC1C6 

15* 

00106 

IS* 

CC106 

17* 

00107 

13* 

CC11C 

19* 

00111 

20* 

00111 

21* 

0C111 

22* 

com 

23* 

00111 

24* 

0011 3 

25* 

30113 

25* 

00113 

2 7* 

00113 

23* 

C0125 

29*' 

00115 

30* 

00115 

31* 

00115 

32* 

CC117 

33* 

00117 

34* 

00117 

35* 

30117 

36* 

00121 

37* 

00121 

33* 

001 21 

39* 

00121 

40* 

“00121 ” 

41* 

00121 

42* 

00121 

43* 

00123 

44* 

00141 

45* 

00142 

46* 

00147 “ 

47*' 

00153 

48* 

D0157 

4 9* 

00166 

50* 

OOl 67 

'51* 

00172’ 

52* 

□ 0173 

53*" 

00174 

54* 

”00174 

50* 


C 

C 

C 

c 

c 

c 

c 


c 

c 


SUCCESSIVE ACCELERATED RTPEAOEMENT METHOD 

NON EQUILIBRIUM FLOW IN THE MERGED SHOCK LAYER OF A BLUNT BODY 
NEAR THC STAGNATION REGION USING FULL N.'S . EQUATIONS 
SLIP CONDITIONS AT THC HALL 
CATALYTIC OR NONCATALYTIC WALL" 


DIMENSION UUOO) fVUOD) tT(lOO) *H(10Q)fPC 100 ) * P2 ( 100 1 » RTI 1 00 ) * 

1D( 100) *RAD (IOC) *X (ICO) tC 02 (100) .CN2I100 )» CN0(100 ).CN(1CC )*C0(100 ) , 
1CN0IU00) tCCL(lOO) *R0(100I 

DIMENSION W0 2(100> *WN2~(100> *NN0(10Q) *WN(100 )»W0( 100 ) * WNOKlfftT) • 

1 WEL ( 1 00 ) » PC L (IOC) 

DIMENSION W2 (ICO I .HO 2 (100) tH 0(100 ) *HN( 100 J *HNO( 100 J.HN0K100) 
DIMENSION HEL(IOO) 

DIMENSION Cl (10) • C2 ( 1 0 ) « C 3 ( 1 D ) * PP1 (10 )• W ( 10 >»CL ( 10 ) 

INPUT HANSEN * S"“V IS C OS I tT^dTT A~ 

DIMENSION XMUSUUOO)tTABS(lOO) «XMUR AT ( 100 ) .X MU ( 100 ) 

DIMENSION TEMP (37) *TAB1(37) *TAB2(37) *TAB3(37)*TAB4(37 ) 

DATA TEMP / 2000 .0 *2 500 . 0. 3000. 0 • 3500 .0 * 4 000 . 0* 4500 .0* 500 0 .0* 

1 5500.0*6000. 0*6500. 0*7000.0*75 00. 0*8000. 0 * 8500 .0* 900C.C*9500 .0 * 

„ 2 10 000.0*10 500.0*11000.0*11500.0*12000.0 *125 00.0*13000.0*13500.0 . 

3 1 4C00.C* 14500.0** 15000.0 *15500. D*1GC0C. 0*16500.0 *17000.0*17500 VO . 

4 13000.0*13500.0*11000.0*19500.0*20000.0/ 

DATA TAB1 /C .886 . 0 .84 6* 0 .830* 0 . 81 5*0 . 80 3* 0 . 792 *C .782 • C'. 773 » 0 .764'* ' 

1 0. 757*0.750*0.743*0.737*0.731*0.725*0 .720*0 .715*0.710*0.706* 

2 0.701*0.6 97*0 .6 93 * 0 .6 89 * 0. 685 *0.6 81 • 0 .67 7. 0.672 .0 .66 8* 0.664 * 

3 0.660*0. 656*0.651*0.647*0.6 43*0 .639 *0 .635*0 .6 30 / 

DATA TAB2 /C . 7 4 2 . 0 . 7 05 • 0 . 6 rTfUTeSTT STD'. 628* 0 . 608 7Er.TTr*"r.'575 • C . 5 6 1 * 

1 0.548*0.536 *0.524*0.514*0.504*0 .495*0.4 86*0 .478*0.470*0 .463* * 

2 0.456*0.448*0 .443*0.437*0.431* C. 426. 0.420* 0 .41 5*0 .407* 0.401 . 

3 0.395*0. 339*0.384*0.376*0.371 *0.365*0 .359*0.352/ 

DATA TAB 3 /C .4 68 . C .457 . 0 .44 5. 07434 *0 .42 3. 0.41Z". 0;T)01 • 0. 1ST* 0.180 * 

1 0.366*0.353.0.342*0.331*0.321*0 .313*0 .304*0 .297*0.290*0.233* 

2 "0 .281**0 .Z7D*'C72E5~tTr.ZET*D^256»0.ZbZ *U. 24 /* U. 243*U. 23 6 ‘ * 0 .23 0 * : 

3 0.224*0.218*0. 213*0.20 6*0.201*0.196*0.1 90*0.185/ 

DATA TAB4 /123 .30* 122 . 20 • 120 .*80 . 1 1 8.10*114 . BIT. 109.70* IZO.OOV" 

1 89.90*75.60*64 .50* 55.70*48.60*42.80*37.90*33.80*30.40*27.40* 

2 24. 9C. 22. 70 *20.80 *19.09. 17. ED*. IE. 27# IS .1 0TT4 .aA7T3YD9. 12 .241 

3 11 .25*10.38*10 .4 8*10.01*9.50* 9.01* 8.75* 8.01*7.51*7.68/ 


READ (5. 70 2) AB1* AB2* AB3* AB4# AB5.AB6. AB7 • A B8 * ABB* AB1 Q* AB11 • AB12 

702 FORM AT ( 1 2A6) ‘ 

READ (5. 700) X BOO* XCA T * XVIS 
RE AD (5 *7 OCT) CASC.CFFR1 
REA0(5*700)T0L*EPSI 

‘READ (5.700) ALT .RTOB iUFS.TFS.TUK 

700 FORMAT ( 6F10 .0 ) 

READ ( 5*701 TROES 7 

701 F0RMAT(E13.8) 

~ETFr=lffri 

ICASE=CASE 


ootmoc 

000000 
■ -QODDOC 
00000(3 

000000 

oooaao_ 

oooooo 
‘ ‘aoouoo 
aooaoo 
ocoooi 
000001 

odonoi 

ooooai_ 

0000 01 

000001 

000001 
000001 ^ _ 
oooooi 
_gooooi__ 

300001 

OOOOOI 

00C0D1 

OOOOOI 

' ‘ ~ o oumji 

OO OOOI . 

OEEUEI 
OOOOOI 
Ocooci 

OOOOOI 

DODOOT 

OOOOOI 

OOOOOI 

OOOOOI 

P0 C 0 C1 

OOOOOI 

OOOOO I 

OOOOOI 

0D0D01 

OOOOOI 

— * onoooi 

oooooi 

DTJOTTOI 

OOOOOI 

D00U22 

000022 

000032 

000041 

00 0 0 50 

00006 2 

000062 

00007 0 

300070 

00007 2 

: 0 00072 


c 



i II iiiiiii in 


DATE 120877 


01 


MAI N 


00174 

56* 

00174 

57* 

00174 

58* 

00175 

59* 

00201 

60* 

00201 

61* 

00201 

62* 

00201 

63* 

002C2 

54* 

00205 

65* 

00207 

56 • 

00210 

67* 

00211 

63* 

CC212 

69* 

00213 

70* 

00214 

71* 

00214 

72* 

DC214 

73* 

00214 

74 • 

00214 

75* 

0021 5 

76* 

00217 

77* 

. 00217 

73* 

00220 

79* 

0D222 

30* 

CC 223 

81* 

00225 

32* 

00226 

S3* 

00230 

34* 

002 30 ' 

85* 

00231 

35 • 

00233 

87* 

00233 

38* 

C0234 

89* 

00236 

90* 

CD 23 7 

91* 

00241 

32* 

00242 

93* 

00244 

94* 

CC244 

95* 

00245 

96 • 

0074 7 

* “97* 

00247 

93* 

00250 

99* 

00252 

100* 

00253 

101* 

00255 

102* 

00235 

103* 

00255 

104* 

00255 

105* 

00256 

106* 

00260 

107* 

00260 

109* 

~C0260' 

“"105*- 

00261 

110* 

30233 

111* 



BOUNDARY CONDITIONS AT THE BODY 

ULTUUTZ 

00 0 07 2 

ROAD ( 5 .7 00 )U(1 ).T(1) 
V€1 >=D - 

000101 

C 00110 

000110 

CL R CPRC$ENTS LEWIS NUMBER AND A IS AVAGACRO NUMBER 
00 2 K-l » 7 
1 CL IK ) -1 • 4 
ELI -1 • 4 

0001115 

000114 

000114 

000116 

CL 2-1 • 4 
P R-Q .75 

0015 117^ 
000120 

GAMMA-1.4 
SCI -PR/CL1 

000122“ 
00012 3 

SC?=PR/EL2 

0013125 

000125 

000125 

PRINT PROGRAM TITLE AND INPUT DATA 
PRINT 10 

uu.Ul/5 

000130 

1C FOPMATtl HI »///25X. ’HYPERSONIC. MCRGED SHOCK LAYER. STAGNATION “LINE" 
2 FLOW FIELD FOR SPHERE OR CIRCULAR CYLINDER’) 

JUU134 

000134 

PRINT 11 

11 FORMAT (/SOX. ’ JA IN-KUMAR-HEND Rlt KS PROGRAM’) 

"" UUU134 

000140 

PRINT 12 

12 FORMAT (45X. ’ (OBTAINED FROM BILL HENDRICKS. JULY 1975)’) 

“[TO 01*40 
000144 

PRINT 13 

13 FORMAT (////40X. ’FEATURES: SIMILARITY SOLUTION OF NAVI! R-STOKES EQU 

000144 
00015 0 

1 AT IONS ’ ) 
PRINT 14 

uUUim; 
00015 0 

14 F0RMATC5DX. ’NONEQU iLIERlUM CHEMISTRY’' 
1 / SOX . ’SLIP BOUNDARY CONDITIONS’) 

HD0I34 
00015 4 

PRINT 19 

19 FORMAT (//31X. ’AVAILABLE OPTIONS:’) 

0110X34 
00016 0 

PRINT 21 

21 FORMAT (50X.35H1 . BODY: (A) SPHERE OR (B> CYLINDER) 

PRINT 22 ' ~~ 

22 FORMAT (SOX. ’2. WALL CATALYTICITY1 (A) NONCATALYTIC WALL OR (B) F 

’ 000160 

00016 4 
000x34 ‘ 

00017 0 

1ULLY CATALYTIC WALL’) " ' ----- - ' 

PRINT 23 

D00170 
00017 0 

X3 FO KM AT (FOX » *3. VlSCOSirYl lAl bUIHLRLAND OR (B f~ HAN SEN ’ ’ S HIGH TE 
IMP. MODEL ( N A CA TR R-50)’) 

000174 

00017 4 

PRINT 24 ‘ 

24 FORMAT (//35X. ’MODIFICATIONS:’! 

000174 

000200 

PRINT 25 

25 FORMAT (50X. *1 . 9/2/76 (KEN JOHNSTON) ENTHALPY COMPUTATION MODIFIED 

000200 

000204 

1 TO ALLOTS FOR TIO? TIALLY 1 */ 55X . 'L XCITED V J.L KA1 lON A~L S 1 A IX T VI NCt NT! 
2 AND KRUGER. P.135) STATEMENTS 30.130.’ /53X. ’PREVIOUS ASSUMPTION 

000204 

000204 

3: FULLY EXCITED VIBRATIONALTTATr’l * ~ ' 

PRINT 26 

00020 4 

26 FORM AT (/ 5DX . *2 • 9/2/76 (kCN JOHNSTON) OPriON ATJUXU T0CTJMPU1E VI3C 
10SITY BY HAN SEN ’ ' S HIGH TE MP.V55X. • MODE L. CSUTHE RLAND ” S FORMULA 

00U21U 
000 210 

/BECOMES - INACCURATE AT HIGH lEHP.I’J 
PRINT 350 

0002 10 
000210 

350 F0RHAT7/5DX.’3. ^)/y/7b (KEN jumNsTON) cAIaLYIIC AND NCNCATALYTTC P 

0002 14 



<! 

Oi 


M.AIN 


DATE 1208 77 


00263 
0C263 

00264 
C02 66 
00266 
C0267 
00305 
CC306 
00310 
CC311 
00313 
C0315 
00316 
C0317 

00321 

00322 
C0324 
0032 6 
00327 
00 330 
00332 
C0333 
00335 
CD337 
00340 
C0341 
00343 
0C344 
00354 
C0354 
00354 
00354 
00354 
00354 

00354 
“ 00354 

00355 
CD363 
00363 

00363 

00364 
003 715' 
00370 
0037C 
00370 
00370 
00370 
00371 — 
00372 
00377 * 

00374 

00375 ' 

00376 

00377 

00400 

00401 “ 


112 * 
113* 
114 * 
115* 
116* 
117* 
113* 
119* 
120 * 
121 * 
122 * 
127* 
124* 
125* 
126* 
127* 
123* 


142* 
143* 
144* 
145* 
143* 
"14 7* 
143* 
149* 
150* 
1 51* 
152* 
‘153* 
154* 
155* 
156* 
157* 
153* 


159* 

160* 

161* 

162* 

TOT* 

164* 

“TD5* 

166* 

167* 


1R0SRAMS C0M8INCD. ’ / 55X » ’(OPTION AVAILABLE TO CHOOSE DESIRED WALL C 

2 AT ALY TIC CONDITION ) * ) - - - - - * 

PRINT 3 51 

351 FORMAT (/50X. *4 . 9/9/76 C KEN JOHNSTON) SPHERE AND CYLINDER PROG RAFTS 
1 COMBINED .’ /55X. MOPTION AVAILABLE TO CHOOSE DESIRED BODY)*) 

PRINT 1 8 t A 81 » A B2 » A C3 i AB4~t AB5 » A B 6 * A B? rAB8"7AB 97 aB T C7X6TI DSBT2 

18 F ORMAT ( ///3 X • ’PROJECT l *»12A6) 

PRINT 27 ~ ' 

27 F0RMAT1/13X. ’OPTIONS USEDIM 
IF CXB0D-0.9 .GT. 0.2) GO TO 29 
PRINT 23 

28 F0RMAT127X. ’SPHERE*) ‘ 

GO TO 41 

29 PRINT 40 

40 FORMAT ( 27 X » ’CIRCULAR CYLINDER* ) 

41 IF (X C AT-C • 9 .GT. 0.2) GO TO 43 
PRINT 42 


000214 
— D00214- 
000 214 
“'000220 
000220 
DO 02 20 
000240 
000240' 
000244 
000244 
000251 


□00235 
000255* 
000257 ” 
00026 3 
PC0263 


129* 

42 

FORMAT (2 7X *’NONCATALYTIC WALL*) 




000273 

1 30 * 


GO TO 45 




00027 3 

131* 

43 

PRINT 44 




D00275 

1-32* 

44 

FORMAT (27X. ’FULLY CATALYTIC WALL’) 




000301 

133* 

45 

IFIXVIS-0.9 .GT. 0.9) GO TO 47 




OD0301 

134* 


PRINT 46 




000305 

135* 

46 

FC RM AT ( 2 7X » ’SUTHERLAND VISCOSITY FORMULA’) 




000311 

136* 


GO TO 49 




000 311 

137* 

47 

PRINT 48 




000313 

139* 

43 

FORMAT t 27 X * ’ HAN SEN * *S HIGH TEMP VISCOSITY FORMULA*) * 




000 317 

139* 

49 

PRINT 15tICASE.tFFRl.um »T(1) .TOL.EPSI 




' ~ 0003TT “ 

140* 

141* 

15 

FORMAT ( //13X. ’RUN N0.:’»I4* 

1 ///15X.* INPUT D ATA'I r ' ' 




000330 

filllM ~ ‘ “ 


1 //18X. ’EFFECTIVE RADIUS (EFFR) = *#F7.4»’ (STANDOFF DISTANCE/BO 

1DY RADIUS)* ' ~ ^ 

2 /3GX.’U1(1) 3 ’ . F3 .5t * (U1 BAR (1 ) /F REEST RE AM SPEED)’ 

3 / 36X * *T 1 ( 1 ) r *.F8.5»* ( T1 BA R fl J/FRETSTREAM STAG. TEMP. T* 

4 /4 X * ’ HEAT TRANSFER COEFF. TOLERANCE (TOL) = ’F8.6. 

5 / 37X * *EPSI - ’ . FST76T7) : 

PRINT 16.ALT.RAD0.UFS.TFS 

16 FORM AT (2 7X. 17 HALT ITU DC (ALT) r ,F9.4.“4H KK/Z3X iZIHBDDY RADIUS “TR - 
1ADB) = * F 9.4 v 4 H C M/1 9 X. 25HF RECST REAM SPEED (UFS) = tF11.3* 8H C 
2M/5EC/13X. 31 HFREESTREAM TEWPERATUKS “nTS7~=— iF?T4 iTZW~ DEG K L LV I N )~ 

PRINT 17.R0FS.TWK 

17 FORM"AT(IOXT*TRrrSTREAH DENS IT Y (RO F S ) ~ = "* »L1 0 .5t ■ — ETI /CH»*3 f 

1 /1 9X» ’ WALL TEMPERATURE (TWK1 ='*»F9.4t’ DEG HELVIN’) 


000330 


000330 
'000330 * 


000 330 
-000330 
000330 
“000340" 
00034 0 
0 0 0 340 
000 340 


COMPUTE AND PRINT FREESTREAM CONDITIONS AND SIMILARITY PARAMETERS 

“V SOUND- TUUDtT* (1 .4 *2 B 7 ,U8*ThS ) » *U . b 

RTFS=(1 .458E-05*TFS**1.5)/(TFS+110.4) 

AM ACH-UFS/VSOUND 


STA GFS- TFS * ( 1 .0*( (G A HMA-1 .0) /2 .0 ) *A MACH* *2 ) 
-RTSrA6=(l.'453E " C5 * ST A GF S»»l.S)/(S T AG FS *11 0 .4) 
RTWALL= (1 .4 58E-05*TWK**1.5)/(TWK+U0.4) 
_ NEYr^UFS^R0FS*TnnJF7RTFS 
REYN=REYF*RTFS/RTSTAG 

REVALL=NEYF*RTFS/R TWALL 


0003 4b 
000346 
000346 
00034 6 
000346 
000346 
-0003 46 
000 357 
' 000 3T1 
00037 4 
" 000 4 0 4 
000416 
000 430 
000435 
00P 441 



MAIN 


DATE 12C877 


<1 


C04C2 

163* 


XLAMDAr (1 .2 5 62* R TFS ) / ( R OFS* { 2. 87C 25 +05 *T FS ) * *0 . 5 ) 

000444 

CC4C3 

169* 


XKf.FSrXL AM CA/R ADO 

0004 6 0 

Q04C4 

170* 


CF Sr R7WALL*TFS/(RTF5*TWK) 

00046 2 

DO 40 5 

1 71 * 


VEARr APACH* (CFS/RCYF)**0.5 

C0C470 

C040C 

172* 


TSTARrC .5 * ( STAG FS+T WK ) 

000500 

CC4C7 

173* 


RTOTARr( 1 . 4 F 8C -0 5* TST AR* * 1 . 5 ) / ( TST A R+ll C . 4 ) 

* 3CC5TJ4 

CC41C 

174* 


COT AR-RTSTAR*TFS/(RTFS*TSTAR) 

0005X5 

CC411 

1 75 * 


XKOGr REYF/(CAMMA*AMACH**2*CSTAR) 

DO 05 22 

C0412 

176* 


PRINT 4 2C 

0C053C 

CC414 

2 77* 


4 CC FC^MAU//* COMPUTED FREES TRCAM DATA:*/) 

DC 05 34 

CC415 

173* 


PRINT 4 01 * V 00 ON 0 » AMA CH * STA GFS 

000534 

CC422 

1 72* 


4 C 1 FCPPAT (7X. f FRCCSTRCAM SPEED OF SOUND (VSCUND) r *tFlC.3,» CM/SCC* 

’QC05“43 ' ‘ _ 

00422 

1 3C* 


1 /11X» 'FREESTRCAM MACH NUM3CR (AMACH) r • ,F6 . 3. /I 2X * • FRCCSTRCAM S 

000543 

CC42Z 

131* 


2TAC T CMP (STAGES) r 'tF?„3» * DCG KELVIN* ) 

000543 

00423 

132* 


PRINT 402 » RTF St RTSTAC 

000543 

CC427 

137* 


4C2 FORMAT {1 4X * 'FRCCSTRCAM VISCOSITY (RTFS) r 'iCIO.Si' CM/CM-SEC* 

0DD551 

00427 

1 34 * 


2 /7X * 'FREEST REAM STAC VISCOSITY (RTSTAG) r ••C10.5» t GM/CM-SECM 

00055 1 

CC42C 

135* 


PRINT 403 

OTJCSBi 

CO 4 32 

135* 


4C3 F0RMAT(/4X» 'SIMILARITY P ARA MET CR S : * / > 

000555 

CC433 

137* 


PRINT 4C4»RCYF*RCYN*RCKALL*XKNFS»V3AR»XKSQ 

DC 05 55 

00443 

133* 


404 FORMAT (3X.* FREE STREAM REYNOLDS NUMBER tREYF) - NFS. 2* 

00056 7 

CC443 

189* 


1 /EX * 'STAGNATION REYNOLDS NUMBER (REYN) r **F9.2* 

DJD 05 67 

0C443 

1 90 * 


2 /12X»*WALL REYNOLDS NUM3ER ( RE WALL J = *»F9.2» 

0C05S7 

CC443 

2 31* 


7 /8X * * FR EESTRC AM KNUDSEN NUMBER (XKNFS) r **F9.4« 

00C5FT 

CC443 

132* 


4 / 37 X * ' VO A R r *,F9.4. 

000567 

CC443 

193* 


5 /25X.’K SQUARED (XKSQ) = *»F9.4»//> 

300567 

00443 

134 • 

- 


000567 

DC443 

195* 

C 

NOTE THE SULST ITJJTION RE Y NrRCYF FOR REMAINDER OF PROGRAM 

300567 " * 

00444 

135* 


REYNrRCYF 

000567 

0C445 

197* 


PRINT 13C ----- - . • • 

0C0571 

00445 

133* 

£ 


000571 

0C445 

199* 

c 

** ** ********************** ********************* * * 

' 000571 

00445 

220* 

c 


000571 

0C44 7 

201* 


UN IRrl 87 

D 00575 

00450 

202* 


ArG .C222*1C.**23 

000577 

C0451 

?C7* 


WTN2-28.C " " ’ ‘ ‘ 

™ 0006 or 

C0452 

204 • 


WTC 2-32 .0 

000603 

0C453 

2Cf * 


WTN=14.C 

D 006 05 

0C454 

225* 


WTOrlG.C 

000607 

CC455 

207* 


WTHC-30. C 

* 3006 IT 

0045G 

208* 


WTNOIr3C.O 

000613 

0C457 

2C9* 


WT CL=1 ./182C. " 

"00 06 14 

00460 

210* 


NDlVrSC 

000616 

CC461 

211* 


DI VN-NDIV 

' "000020 

00452 

212* 


DCLTArl ,/DIVN 

000623 

00463 

21 3* 


NDIV1-NC IV+2 

000626“ 

00464 

214* 


HAL Tr25 *0 

000633 

CC464 

215* 

c 


” 000633 

00464 

215* 

c 


□0063 3 

00464 

217* 

c 


■00 0633 - 

00465 

213* 


READ 30 C» (U (I ) r V (I) »T(I) »RO(I) » P2 ( I ) »I=2*NDIV1 

000635 

00477 

219* 


READ 3DO.(C02(I).CNOI(I)*CNOlI)tC7imiCO(Tl.T=2*NDXV) 

000657" " 

004 77 

220* 

c 


00065 3 

C0477 

~221 * 

c 

************** ******************** ififnrrriTffi * 

0006 53 

00477 

222* 

c 


000653 

00477 

. 223* 

c 

BOUNDARY CONDITIONS AT FR ECSTRLATT “ " “ 

000653 - 



<1 

00 


MAIN 


DATE 12C677 


0C511 

224* 



U ( N DIV1 )= 1 • 

00067 0 

CC 512 

225* 



VC NDIV1)=-1 . 

3CC672 

00 513 

226* 



TCN0IV1 )=TFS/STAGFS 

00067 3 

CO 51 4 

227* 



RC CNDIV1 1 = 1 . 

Q 006 76 

00515 

228* 



RT ( NDIV1I=1 .0 

000677 

CC51& 

229* 



P2CNDIV1 >=0. ' "" ' 

CTTCTOC 

00517 

230* 



CN2 INDIV1 l=C .767 

000701 

0C52C 

231* 



C02(NDIV1)=C.233 

0007D3 

00521 

232 * 



CNC (NO I VI ) = 0. 

000705 

CC 522 

23 3* 



CO (NDIV1 ) = C. 

000706 

00523 

234 • 



CNC NDIV1 1=0 . 

000707 

CC524 

235* 



CNCICNDIV1 »rC. • ' ' 

“OBfiVTC * 

00525 

23 5 * 



CCLCNDIV1 ) = C. 

000711 

CC525 

237* 

C 



00C711 

OC 525 

233* 

C 



000711 

C 0 52 5 

239* 

c 



000711 

00525 

24 0* 

c 


CATALYTIC BOUNDARY CONDITIONS 

000711 

CC 525 

241* 

c 


_ 

000711 " 

0C525 

242* 

c 



000711 

CC526 

247* 



IFIXCAT-C.9 .LT. 0.2)G0 TO 450 

0 0 07 12 

C053C 

244 • 



CN2 (1 >=0.757 

0C0717 

C 0 53 1 

245* 



CO 2 ( 1 >=0.233 

00072C 

CC 5 32 

245* 



CN0C1 >=C. 

000721 

00533 

247* 



CN C 1 >=0.0 " 


0053*1 

248* 



C0( 1 > = C.0 

00072 3 

CO 535 

249* 



CNCI C 1 )- CNOI (2 ) 

000724 

00536 

250* 



can ) = cnoi( 2 >*wtel/wtnoi 

000726 

CO 53 7 

251* 



PELCl >=AtRO(l)*ROFS*CNOICl)/WTNOI 

000732 

CC54C 

rnran 

252* 

*>C7* 

n 


GO TO 451 

000740 

L L jHU 


L 



' IJTTOTW" 

00 540 

254* 

c 



000740 

CC 54 0 

255* 

c 


NOKCATALYTIC WALL BOUNDARY CONDITIONS (FOR FIRST ITERATION) 

000740 

00540 

255* 

c 



000740 

CO 54 0 

257* 

c 


***************************************************** 

0007 40 

00541 

253* 


450 

CN2 C 2 ) =1 .0-C02(2)-CN0 C2 )-C N( 2>-C 0 (2 )-C NOI ( 2 1 • (1 .Q«U TEL / WTNOI ) 

0007 4 2 

CC 542 

259* 



CN2C1UCN2C2) ~ ' “ — 

0CC754 

00543 ' 

260* 



C 02 (1 )=C02(2) 

000755 

0 C 54 4 

261* 



CNC ( 1 ) =CNO (2) 

000757' 

00545 

262* 



CNC 1 ) = C N ( 2 ) 

000761 

OC 54 6 

267* 



CO C 1 ) =C0 ( 2 > 

‘ CC076T 

CC547 
nor cp 

264* 



CNOI Cl ) =CNQ I C 2) 

p ri f i i —aii A T i ff % f fTn - r i rriTffT — ■ ■ ■■■ ■ -■»■■■- — »■< — - 

000765 

UL DDL 

£ bb * 



CLLl 1 J-CN0II2 ) *KTlL/wTNOI 

DTTT7F7 

00551 

266 * 



PEL Cl >=0. 

00077 3 

CC 551 

267* 

c 



' 000773" 

00551 

268* 

c 



00077 3 

CC551 

269* 

c 


- * 

000773 

00551 

270* 

c 

451 

INITIAL GUCSS FOR CNOI*CEL A NO H 

00077 3 

CC 552 

271* 


KIN2-C .0 

000775 

00553 

272* 



HI02=C. 

000775 

CC554 

273* 



HIN0=305C.*1C.**7 ' - - - 

'0 0 07 76 

00555 

274* 



HIN=33620.*10.**7 

aoiooo 

' 00556 

275* 



HIC=1543C.*10.**7 " ' ' - • - " " ‘ “ 

001002 

00557 

276* 



H IN 01=9000. *10. **7 

001004 

~DO 56C 

277* 



H I CL — Co 

'"001006 

00561 

2^8* 



W1=8.3143*1C.**7*STAGFS/UFS**2 

001007 

CC562 

279* 



DO 30 1=1 tNDIVl -- 

001200 



PAIN 


DATE 12 0877 


00555 

230* 

CN7C I) -1 .-C0 2 (I)-CNO II) -C N II)- C 0 (I) -C NOI Cl Mll.+U TEL/4 TNOU 

001216 

0C566 

281* 

CELCD-CNOKIJ *WTEL/WTN0I 

nor ? 26 - 

00567 

233* 

W31 =T 1 1 )*STAGFS 

001232 

CP57C 

287* 

VN 2=3390 ./W31 

001235 

00571 

294* 

V 02 =2270* /X 31 

001240 

C0572 

235* 

VNC = 2 740 ,/W31 

CC1243 

09573 

235* 

VN0I=274Q */V!31 

0C124 6 

CC574 

287* 

IF (VN 2 .GT. 8C.OGC TO 455 

00 12 47 

00576 

233* 

VI3=CCN2< 11/28. )*VN2/(EXP(VN2)-1 • ) + 1 C02I 11/32. )*V02/CXP( V02 1-1 • 1 

00125 2 

0C576 

789* 

1* t CNO (I>/30.)*VNO/ IEXP(VN0)-1.)+<CN0I(I 1/30* 1 *VNOI /( EXP ( VNOI 1 -1 • ) 

901252 

C0577 

290 * 

GO TO 456 

CC1317 

C06CC 

291 * 

455 VILrO.C 

901321 

00601 

292* 

455 CONTINUE 

001322 

CD6C2 

297* 

H( I)=3.5*(CN2 1 D/26.+C02 (I)/32*+CN0(I J/3C.+CN0II I ) / 3C# ) + VI6 

001322 

00602 

2 04 * 

142.5* (CNII)/14.+C0t I)/16.*CN0I(I)/3C. ) 

001322 

CC6C3 

295* 

35 H ( I ) = W 1 * T I I > * H ( I ) 4 ( C NO (I)*H1N0+CN(I)*HIN+C0(I)*HI0+CN01(I)«HIK'CI) 

001347 

00603 

293 • 

1 /ur 5**2 

001347 

00603 

297* 

C 

— ' DTIT 47 

crcr3 

799* 

C 

001347 

0C6C5 

390* 

DO 5GC Kr 1 ,2000 

00137 4 

C'-EIC 

301* 

Q= EF FR *0 EL T A 

901377 

0C610 

302* 

C TO FIND P AMD HUE AT ALL POINTS 

00 1 37 7 

00611 

3n* 

CO 50 1= 1 * NC IV1 

-joivcEr 

006m 

324* 

T AD S (I ) = T 1 1 ) * 57A GF5 

001406 

crcis 

705* 

IF (I .EG. 1) GO TO 31 

901411 

0061 7 

305* 

P( I >=R0 (I)*T (II*W1* I CN2 (IJ/UTN2 + C02 ( I )/WT 02+ C N0( I ) /M TNO+ CN( I ) /M TN 

0014 14 

00617 

3r7* 

l+Ci/II)/WT0+2.*CN0I (IJ/WTNOI) 

901414 

00620 

303* 

31 XMUSU(I) = 1 .453C-05*TA3S* I)**1.5/« TAB5C II + 110.4) 

0014 4 4 

C0621 

309* 

IFIXVIS-C.9 .LT. 0.2)00 TO 33 

0014 57" 

00623 

310* 

IFITA3SU) .GE. 2000.0) GO TO 32 

001461 

C0625 

711* 

33 XMUR AT (1 1 = 1.0 

0C1466 

0C626 

312* 

GO TO 37 

001467 

00627 

313* 

32 T A lS =TA6 S 1 1 ) 

901471 

00630 

314* 

CALL INTERP(1»4»1.37*TEMP.TABS) 

0C147 2 

CC631 

315* 

CALL INTERP(2*4*li37»TABl»SN2N) 

~0Cl5"02~ - " ~ 

00532 

316* 

CALL INTERP I 2 t4 »1 »37*TA92»SNNI 

001512 

00 633 

717* 

CALL INTERPt2.4.1.37rTAB3.SNE) 

901522 

00634 

318* 

CALL INTERP(2.4.1.37iTAB4fSEEJ 

001532 

CD 63 5 

319* 

XLKN2=1 .03329* CN 2 1 1) +.87542*C 02 tl) *2.5310 *CN(I )*SN2N 

301542 

00635 

320* 

1 *2 .120 36 *C0 1 1) *SN2N + .94819*CN0tI) + .94 819*CN0I(I) 

001542 

CCC36 

321* 

XLPC2 = l.C696*CN2(T)4.9T341T*rcr2TrT ; »2".648 3'*CrRmtSN2N ‘ ” 

1701556 

00636 

322* 

1 42.2146*C0(I)*SN2N4.9303*CN0(I)+.9803*CN0I(I) 

0C15S6 

C0637 

72 3* 

XL MN = • 8949* CN2 (I ) • SN2 N + . 7665* C 02 (I ) *S N2 N+ 2 . 0666 *CN 1 1 )*5NN 

0016 12 

00637 

324* 

1 +1 .7508* CO (II* SNN+ .825 9*C NO ( 1 1 * S N2 N+. 8259 *C NOI 1 I ) * SN2 N 

001612 

006*10 

325* 

XLMC=.9159*CN2 ( I > * SN2 N+. 7830* C 02 (I > *S N2 N+Z. 1391 *CN(I )*SNN 

001642 

006*10 

325* 

1 +1 .80 83* CO (I)*SNN*.8444*CN0(II*SN2N4.8444 *CNOI(I ) *SN2N 

001642 

00641 

32 7* 

XLPNO=l.C516*CN2(ri4.3gg9*C'02'fir42'^'gn6*CNrT)*SN2N 

"0 0 16 T 2 

006*11 

328* 

1 *2 .1680* CO CI)*SN2N+.9644*CN0II) +.9644 *CNOI(I) 

00167 2 

006*12 

329* 

XLPNO 1=1 •C516*CN2(I)4.8759*C02(I) + 3.2475*CNCI) *SN2N 

001716 

006*12 

330* 

1 *7 .16 80* CO <I)*SN2N4.9644*CN0<I) + .9644 *CNOI(I) 

001716 

0064 3 

331* 

XLPEL = .7307*CN2 Cl) *S NN+. 6393* C02 I I ) *S NN+l .4 613 *CNCI) • SNE 

• -DC17 31 

00643 

332* 

1 4.2736*C0( I)*SNE + .6819*CN0II)*SNN+.6819*CN0ia)*SEE 

0017 31 

00643 

333* 

2 + 52656. 24*CEL (I)*SCT ' ' - - 

- 0P1T31 

00644 

334* 

XMUR AT (I)=l .016 5*CN2CI)/XLHN24.9509*C02( D/XLM02 

001765 

00644 

335* 

1 +1.4376*CN tI)/XLMN+1.8083*COfI)/XLHO4.982O*CNO(I)/XLKNO 

001765 - 


00 

o 
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0CS44 

3 35 * 


2 ♦ 0 . 98 20 * CNO I( I) /XL MN 01 + 229 . 46 9*CCL (I)/ XLMEL 

001765 

00645 

337* 


37 XMU(I) - XMURATCIJ *XKU5U(I) 

002022 

00646 

333* 


5C RTI I) 3 XMUm/RTFS 

0C2024 

0064 6 

230* 

C 

TO FIND PRESSURE CN THE SURFACE USING MOMENTUM EQUATION 

0C2024 

C0G5C 

340* 


U10=<4. *U(2 >-U(3)-3.*Um)/DCLTA/2. 

002030 

CO 6 5 1 

341* 


V1C=C4.*V(2)-V t31-3.*VCl 1 ) /CELT A/ 2. 

002041 ' 

00652 

342* 


RT 1 C- ( 4 • *RT (2 )- RT ( 3 ) - 3 . *RT (1 ) ) /DELTA/2 • 

00205 1 

0C6E3 

242* 


V2C=-V{4 )+4.*VI3)-5,*V(2)42.*V(l) 

002062 

00654 

344 * 


DEL 4 V- V (5 1-4 . *V (4 ) ♦ S . * V ( 3 ) -4 . * V ( 2) + V(l ) 

00207 4 

CC655 

345* 


V2 C-V 2C+ 11 ./12 .*DCL4V 

002105 

CCC56 

345 • 


V27-V2C/DELTA**2 

002110 

CCG56 

747* 

c 

SPHCRC/CYLINDER OPTION FOLLOWS 

DC 21 1C 

0CG57 

34 3* 


I F( X50D-0 . 9 .GT. 0.2) 30 TO 51 

002112 

CT 65 7 

344* 

c 

SPHERE 

0C2112 

CCS61 

3 30 * 


P12-V2C*RT (1 )/E FFR + VIS * 12 . *R T ( 1 > + RT 1 0/ EFF R ) + »J 1 0 *R T ( 1 > /2 . -Utl ) * ( 

002114 

CC6G1 

351* 


1RT1C+3.5*RT(1) *EFFR) 

0C2114 

00662 

332* 


GO TO 52 

002137 

0C6C2 

357* 

c 

CYLINDER 

002137 

OCSC 3 

354 • 


51 P1C=V2C*RT(1 )/EFFR + V10* (R T ( 1 ) + RT 1 0/ EFF R ) +U 1C *RT ( 1 ) /4 .C -U ( 1 ) • 

002141 

CCS63 

355* 


1 ( FTlC/2 .C + 7.C/4 .C*RT (1) *EFFR) 

0 C 21 41 

C0S64 

355* 


52 P12-P10*4 .C/3.C/RCYN 

002166 

CC685 

357* 


P( 1 ) = P C 2 )-6ELTA*MC 

002172 

0C666 

35 3* 


R0( 1 >=P (1 )/ Wl/T (1 )/ CCN2 (1>/WTN2*C02 (1 )/WT 02+CN0U >/V TN0+CNC1 ) /WTN 

00217 5 

CC 666 

359* 


1 + CCC 1 ) /WTO + 2.* CNO I Cl )/WTNOI) 

OC’2175 “* 

C06S7 

360 * 


DO GO N -2 *N D I V 

0C2252 

CC672 

361* 


EN=N 

3022 67 

00673 

362* 

52 

X I N )=1 • + C EN- 1 .) *Q 

002 27 2 

CO 673 

3G7* 

C 


0C22 72 

0C673 

364* 

c 

SLIP CONDITIONS 

0C2 27 2 

CC 673 

7G5 • 

c 


DC2272 

00675 

356* 


AC=1.C 

0C23C1 

CCG76 

367* 


WC 1)=WTN2 

0023 03 

C0G77 

353* 


W (2 )rWT02 

002305 

CC7CC 

369* 


K ( 7 ) - W T N 

0023C7 

007C1 

370* 


W ( 4 ) -WTO 

002311 

CC7C2 

771* 


KCn=KTNO 

- •- ‘ - PC2313 

00703 

377* 


W ( 5 )=WTNOI 

002315 

CC704 

377* 


W< ‘’J-WTEL 

0C231T 

0C7C5 

374* 


Cl f 1 ) -CN2 (1 ) 

002321 

CC 70 6 

375* 


Cl (2 ) -C02 ( 1 ) 

002323 

00707 

375* 


Cl C 3) = CN( 1 ) 

002325 

0071 0 

377* 


Cl ( 4 ) -CO ( 1 > — 

" " CC2327 

00711 

373* 


C 1 f 5 ) - C NO ( 1 ) 

002331 

GC712 

379* 


Cl C 6 ) -CN 01 C 1 ) 

DC2373 

00713 

3 8C * 


Cl (7)-CCL (1 ) 

002335 

0C714 

381* 


C2 (1 ) -CN 2(2) 

0073 37 " 

00715 

332* 


C2( 2 J-C02 (2 ) 

002341 

* CC71G 

383* 


C2 ( 3 ) -CN ( 2 ) ‘ " ' ’ — 

orara 

C0717 

334 • 


C2 ( 4 ) - C 0 ( 2 ) 

002345 

CC72C 

385* 


C2 1 5 ) -CNO ( 2 ) 

007747 

QC721 

335* 


C2(6)=CN0IC2) 

00235 1 

C0722 

787* 


C2l7jrCEL<2> 

302753 

00723 

333 * 


C 3 < 1 )-CN2 ( 3 ) 

002355 

" C"0 72 4 

38 9* 


C3 (2 )-C02C 3) --- 

307757 

00725 

390* 


C3( 3)=CN(3) 

002 361 

00726 

391* 


C3C4)=C0(3I ... — 

002363 
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0C727 

33? • 


C 3 < 5 )- C NO (3 ) 

00*2 36 5 

CC73C 

333* 


C3 1 C )=CNOI (3) 

302367 

3C731 

334 • 


C 3 ( 7 >-CCL (3 ) 

002371 

02732 

7 or * 


TW -TUK/S TA GTS 

DC2373 

DC733 

3 35* 


SUM 1 -0 . 0 

002 37 5 

CC734 

397* 


DO 61 J- 1 » 7 

DC 24 C6 

30737 

39 3 * 


61 SUMl-SUMl +C1 ( J) /U ( J ) 

OC 24 06 

CC 7*4 1 

399* 


LO 62 J-1.7 

3024 15 

22744 

432* 


C2 PP1 1 J 1 — (Cl ( J ) /W f J) J *P ( 1 ) /SUMl 

DC2415 

30746 

4 p 2 * 


sur 2 =c.c 

332422 

22747 

42'’ * 


DC 53 J -1 1 7 

CC2426 

2272? 

42 


f 3 St‘P3 = SUK2+ (C2 1 Jl-Cll J) >/Q 

2024 26 

32754 

424 • 


6 UP 3-1 . 

C02433 

CO 755 

4C5* 


SU M 4 -C . 2 

302435 

22756 

4 32 * 


DO 35 J=1 ♦ 7 

002441 

CC7S1 

437* 


C r 'Sur 4 rcir 4+ <TW*RC ( 1 )*C1 1 J) / (W t J) *SUM1 ) **1 . 5 ) *( l.+RTf 1 1 /Z. /PP1 ( J >/ 

302442 

23721 

43] * 


! R CYN* ( U fl )-4.*V(2)/Q+V(3)/QM( (3.1416* GA MM A *W(J)*SUMl*TFS/2./T(ll/ 

3C24 4 2 

CC 761 

425* 


2ST AGFC )* *0 .5 )/Cl f J )* ( 2.-ACI /AC* AMACH/RCYN *EL ( J l/PR*( K( J )*SUP1 > • 

002442 

3C761 

412* 


3RTI 2 )/R0( 1 ) * S UM 2 ) 

0C2442 

CC763 

411* 


SUMO-C.C 

302514 

23764 

412* 


DC 66 J-l #7 

0C2520 

CC 767 

41":* 


6 6 CL K5=SUP 5+ Cl ( J ) *R0 <1 ) /( C W C J) *SUM1 ) **1 .5 ) • ( 1 . + RT 1 11/2. /FPK JI/RCYN 

002520 

2C767 

414* 


1 * ( J ( 1 ) - 4 . *\M2 )/Q + Vt 3)/3) > 

C0252 0 

CC 771 

41 5 * 


USLIP = PT (1 >/P t D/RCYN/AMACH* ( I6.2832*STAGFS *T( 1 ) /C AHP A/TFS ) * *T . 51* 

' '332544 

3C7*»1 

413* 


l(2.-AC)/AC/2.*(4.*b(2)-3.*U(l)-U<3) 1/2./Q 

002544 

CC 772 

41 7* 


TSLIF = 12 .-AO/AC* (3.1 436**C.5) 

3C26 00 

C0773 

419* 


TSL IP = T CL IP • (AM ACH*R7.C1 ) •GAMMA/’R/RCYN/? ./ (GAHMA-1 . )*( (GAMMA *TFS/ 

0026 0 4 

CC773 

419* 


12. /T (D/STACFS )** f 5)* (T(2»-T(1J )/0-2.E*Ptll *RT(1 )♦ 

3C26 C4 

□ 0773 

422* 


♦am; ch**3. 

0C26C 4 

CC773 

423* 


2*( (CAMMA*TFC/CTACFS»**1.5)/R0ll)/RCYN/( (2.*T(1 ))**D.5 1*EL1/PR* "" * 

3026 04 

0C773 

422* 


3 3 UM2 ) * SUM 3 

0C2604 

CC774 

423* 


TSLIP = (TSLIP*SL'M4 >/SUK5 

0C2661 

02775 

424* 


USCONT- (USLIP-U(l) 1/U(1 ) 

00266 4 

0C776 

425* 


US CONS-USCCNT 

3C267C 

00777 

425* 


USD ONT = A3S ( USCONT ) 

0C2G71 

01CDC 

427* 


IF(USCCNT.GT.CPSI) U 1 1 ) = U (1 ) +USC ONS/USC ONT *CPS I *U 1 1 ) 

3325 73 

21002 

423* 


IF( USCGNT.LC.CPCI) UtllrUSLIP 

0C2705 

Cl 004 

429* 


TS CON T“ ( TSL IP -T 11 J )/T (1) 

3027X3 

D1CC5 

430* 


T SC ONS- TS CO N T 

0C2717 

C1CC6 

431* 


TSCONT-ABS(TSCONT) 

3027ZC 

0 1 007 

432* 


IF( TSCONT.GT.CPSI) T(1 )-T( 11 +TSC0NS/T5C0NT *C PS I *T ( 1 ) 

0C2722 

01011 

433* 


IFfTSCONT.LC.CPSI> Tm = T$LXP 

3C2734 — 

0101 3 

434* 


RTf 1 )= (1.45 SC-05* (T 1 1 > * STA SFS) *• 1 . 5 ) / t IT ( 1 ) *ST AGF S*ll 0 .4 )*RTFS> 

002742 

C1C13 

435* 

C 

CORRCCTIGN TO V PROFILC 

0C2T4 Z 

01014 

435* 


DO 70 N“2 »N DI V 

00276 2 

01C17 

437* 


V2=V (N*l )-2.*V tN MVIN-l) 

002T64T 

01 020 

439* 


P1-P1N+1 )-P (N-l ) 

00277 1 

01C21 

439* 


VI =V(N*1 l-V(N-l) 

DC2774 

01022 

440* 


RT1 =RT(N*1 )-RT(N-l) 

002777 

01023 

441* 


U 1 -U ( N *1 ) -U ( N - 1 ) 

003302 - 

01023 

442* 

C 

SPHCRC/CYLINDCR OPTION FOLLOWS 

00 300 2 

C1024 

443* 


IF ( X BOD-0 .9 .GT. 0.2) GO TO 67 

003005* 

01024 

444* 

c 

SPHCRC 

003005 

01026 

445* 


CAPV--RT (N )*V2/Q**Z+0.3T5*RTYW(J* ( PI + RO TN 1 *V n'Tr*VTT^VT7G " 

— 003007 

01026 

446* 


1 *( RT(N)/X(Nl+RTl/4./Q)-RT(N)*Ul/4./Q/X<N) 

003007 

01027 

447* 


CAPV = CAPV+ CU(N )4VrFm/Xm*r3.5*RTTN7/XTin+RTl/2./QI ■ - 

• 003035 - 


00 



00 

to 
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01030 

443 ♦ 

DCRV-2. *RT(N )/Q**2+0.375*RCYN/Q*R0( N)* VI +3.5 *RT t N ) /X ( N 1 * *2 

00 3 05 2 

C1C3C 

449* 

l+RTl/2./G/X(N> 

- 00 3052 

C 1 3 31 

4 50 * 

GO TO 68 

00307 2 

01C31 

451* 

C CYLINDER 

003072 

C 1C 32 

452* 

67 CAPV--RTtN)*V 2/ 0**2+0.37 5* REYN/ 3* ( PI + R0( N ) *V t N ) *V 1 J-Vl/Q/2.0 

00 307 4 

C1C32 

45T* 

1 • (RT (N) /X (N )+RTl/2.0/G) -RT (N)*U1/8.C/Q/X (N ) ----- 

003074 

0103 3 

454 * 

CAPV=CAPV+(UCN)+V<N 1) /X « Nl • C7.0* RTC N)/XI NJ +RT1 /0> /4.0 

003126 

01034 

455* 

DERV-2«C*RT (N)/Q**2+C.375*RCYN/Q*R0{N ) *V1 +( 7.0 *RT<N ) /XC N)**2 

• DC 31 45 

CIO 34 

455* 

1 +RT1/Q/XIN ) )/4 .0 

003145 

01034 

457* 

C TO FIND CCCV 

003145 

01035 

453* 

68 SI3 V-A3S(CPSI*DCRV*V(M1 /CAPV1 

0C316 7 

C1C3G 

459* 

IFtSIGV.LT .1 • ) CCCV-SIGV ' “ 

"007174 

01040 

465* 

IFCSICV.GE.l • ) CCC V~ 1 • 

003202 

C1C42 

461* 

?r V(N)-V(N ) -CCCV *CAP V/DCRV 

003210 

Cl 044 

432* 

VC2 lz(V (3 )-4.*DCLTA*UCll*EFFR>/4. 

003216 

01C44 

463* 

C TO FIND CORRECTION TO U PROFILE 

0032 16 

01045 

4G4 * 

DO ac N=2*NDIV 

003224 

C1C5C 

4 G5 * 

U2=U(N+1 )-2.*U (N)+U(N-1) 

' 00 32 45 

□ 1051 

4 55 * 

U1 - U ( N + 1 ) -U ( N— 1 ) 

0C325 2 

01052 

467* 

RT1 = RT(N + 1 ) -RT (N-l > 

0032 55 

C1C53 

4G 3 * 

V1"V(N + 1)-V(N-1 ) 

0 0 3 26 0 

Cl 0 54 

4G3* 

, CAPUz-RT (N >*U2/G**2+RCYN*RC IN)* VlN)*Ul/2./Q 

003263 

C1C54 

475* 

C SPHERE/CYLINDER OPTION FOLLOWS 

0 0 3 26 3 

C1C55 

471* 

IF ( X 00 D-0 • 3 .GT. C.2) GC TO 71 

"DO 7274 

01055 

472* 

C SPMCRC 

00327 4 

C1C57 

473* 

CAPU=CAPU+ CU <N >+V(N> > /X (N)* (8 .*RT< N)/3./X<N J+RT1/2. /Q+REYN*R0(NI 

003276 

C1C57 

474 • 

1 *U< Mil 

003276 

C1CGC 

4 75* 

DEPU=2.*RT (N)/Q**2+8 .*RT <N)/3./X{N>**2+RTl/2./Q/X(N >+FCYN*ROCN > 

903314 " 

01050 

473 * 

1*(2 • *U ( N ) +V(N))/X(N) 

003314 

01061 

477* 

GO TO 73 " v 

‘ ‘ 90 33 35 — 

010G1 

473* 

C CYLINDER 

003335 

C1CG2 

479* 

71 CAFU = CAPU+ (U (N l+VCN) )/XCN>* (7.*RTC N|/3./XCN>+PTl/2. /Q+PCYN *R0 I N J 

307337 

C1CG? 

4 35* 

1 *U(N » > 

003337 

C1CS3 

4 01* 

DCRU=2.*RT 1N)/G**2+7.C*RT IN)/3./X(NJ * *2 +RT 1 /2 . /Q /X (N ) +RCYN *R0 ( N ) 

DC 3362 “ 

010G3 

432* 

1 *(2.0*U(N)+V(N))/X(N) 

00336 2 

01CG4 

483* 

73 CAFUzCAPU+RT(N )*Vl/G.C/Q/TrRT-TJT7T3*rRTT7n^NT+RTT74*;r707— 

0034 04 

01064 

4 34 * 

1 +2.C*P2(N)*REYN/X(N) 

003404 

C10G4 

435* 

C TO FIND CCCU 

-1303464- 

Cl 0 65 

436 * 

SICUzABS(.CPSI*DERU*U(N)/CAPUI 

0034 31 

C1CGG 

487* 

IFCSICU.LT .l.)CCCU=SIGU 

" 603437' 

01070 

433* 

IFISIGU.GE.l.)CCCUrl. 

003445 

CIO 72 

489* 

3 C UC N) zl! IN ) — CCCU *C AP U/DCRU 

003453 

01072 

4 90 * 

C SPECIES PRODUCTION RATE 

00345 3 

C1C74 

491* 

DO 90 Nz 1 » NO IV 1 

9074 01 

C1077 

432* 

U5zTtN>*STAGFS 

00347 4 

CHOC 

493* 

WG ZUNIR* W5 

“097476 

01101 

494 * 

W7ZH 8000 ./WG 

003500 

Cl 102 

495* 

W8ZW7**1.5 

0035 03 

01103 

493* 

W9Z224900./WG * 

003510 

01104 

497* 

W10Z15OCOO./W6 " 

00 3513 

01105 

493 * 

WllzuiC**2 

003516 

Q31C6 

499* 

W12Z39100./WG “ ‘ ' 

0035 20 — 

01107 

500* 

W13Z7C80./W6 

003523 

'01110 

501* 

W1 4=75500 • ?N6 

00352b 

01111 

502* 

W15Z128500./W6 

003531 

"01112 

5C3* 

W1GZ85520./WG *“ 

093534 
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01113 

504 * 


U17-6329C ,/VS 

003537 

C 1 1 1 4 

505* 


U15=W5**0.£ 

0C35 42 

01115 

505* 


W1 7 = W5 **2 .5 

003547 

Cl 11 6 

507* 


W2C=5CC./W6 

003554 

C 1 1 1 7 

503* 


U2:0=-1C./W6 

0C3557 

ci: 2 r 

509* 


IF (V.’ 7 • C7 .8C. )C0 TO 61 

003562' 

31122 

51 C * 


FK1=2*5*10. **11*W13*W3/EX?CW7) 

003566 

C 1 1 22 

511* 


GO T 0 82 

CO 35 76 

Cl 1 24 

512* 

31 

FK1 =C • 

0G360Q 

C1125 

51 7 * 

32 

RK1=2.CP*1C.** C*W5*K3 

QC3G 01 

C112S 

514* 


IF( W3.GT.3C. JGO TO 33 

CC36C4 

C2l3r 

515* 


FK 2=1 . 7* 10 •** 1 2* W1 B* W 9/CXP 1 WS ) 

0 D 36 1C 

01131 

5 1 C * 


GO TO 34 

003620 

C1132 

517* 

C 7 

FK2=C. 

003622 

Cll 33 

513* 

34 

RK2 = 9.444 *1 C.**10*W13*W9 

003627 

Cl 1 34 

519* 


IFIWir.CT.3C.] GO TO 85 

CD 36 26 

Cl 1 36 

5 ? C * 


FK7=7.*lC.**ia*W13*Wll/EXP<W10 * 

0C36 3 2 

Cl 1 37 

521* 


CO TO SB 

003642 

Cl 1 40 

527 * 

35 

r K 7 =0 • 

003644 

Cl 1 41 

523 * 

8G 

RK3=1.75*1C.**1C*W18*K11 

CD 36 45 

Cl 1 42 

574* 


IF* W12.CT.8C. )G0 TO 87 

C0365 C 

Cl 1 4 4 

525* 


FK4 = 3.2*1C .**9*W5/CXP <W12) 

3C3G54 

Cl 1 45 

52 0 * 


GO TO S3 

QC3663 

C114G 

527* 

07 

FK 4 = C • 

003665 

Cl 147 

573* 


IFC W13.3T.8C) GO TO 391 

003665 

Cll 51 

529* 

33 

RK 4=13.33*1C.**2*WS/CXP(W13) 

0036 72 

Cll 52 

530* 


GO TO 932 . 

CD370D 

Cll 53 

571* 

831 

RK 4 = C • C 

0037 02 

01154 

532* 

332 

IFf W14 .CT .30. JGO TO 91 

OC3703 

C115G 

533* 


FK r = 7.*lC.**13/CXP (Ml 41 

0C37DS 

C 1 1 57 

534* 


GO TO 92 

0C3714 

Cll EC 

535* 

a i 

FK "=C. 

CC3716 

31161 

5 36* 


IF* W2C.CT .30. ) GO TO 921 

0C3716 

Cl 1 63 

577* 

32 

RK r - = l .55*1C.**13/CXP * W 2 C ) 

3C37 2? 

01 1 64 

53 3* 


GO TO 922 

0C373C 

C11C5 

539* 

521 

RK f = C . G 

DO 37 32 

Cl 1 55 

540* 

322 

IF* W15.GT .30 . ICO TO 93 

0C3733 

Cl 170 

54 1 * 


FKr = S.l* 10 .**2 4/W19/[IXPlW15) 

0C 37 36 

Cl 1 71 

542* 


GO TO 34 

003745 

Cll 72 

543* 

93 

FK C=C • 

003747 

01173 

54 4 • 

94 

IF* W16.GT.8C.)G0 TO 95 

00375 0 

Cl 1 75 

545* 


RK6=4.79*1C.**23/W19/CXP ( W1 6 ) — ' " 

. CTC7T5T- 

C1176 

546* 


GO TO 96 

00376 2 

C 2 1 77 

54 7* 

9 r 

RK C = C • 

OD 3764 

C12CG 

54 3* 

36 

IF* W17.CT .30. JGO TO 97 

C03765 

C12C2 

549* 


FK 7=6.4* 1C. **9*K18/CXPCW17) 

003770 

012D3 

55C* 


GO TO 93 

003777 

C12C4 

551* 

97 

FK 7=C • ' ..... 

• - ~ ’ ’004001' 

Cl 20 5 

552* 

93 

RK 7 =1 .73*10 .**1 9/W5/CXP (W210) . 

004 00 2 

C12C6 

553* 


W2=R0*N) *ROFS 

0C4D1O 

01207 

554 • 


W3= W2* *2 

004 013 

0121C 

555* 


R1 =W3* ( CN <N J/WTN +CNO * N J /UTNO + 2. *CN2 * N )/ WT N2+9. *C 02 (N ) AT02 + 

trc4tri5 

01210 

555* 


125.*C0(NJ/WT0) 
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